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The  surface  waves  from  mine  blasts  recorded  on  the  Central  Minnesota  Seismic  Array 
were  analyzed  for  propagation  characteristics  and  velocities.  The  seismic  array  was  a  six 
station  vertical  array  with  a  diameter  of  28  kilometers.  It  was  located  in  east  central 
Minnesota  and  was  emplaced  in  1976  to  monitor  local  seismicity.  It  routinely  recorded 
mine  blasts  which  occurred  in  the  Mesabi  Range  of  northern  Muuiesota,  about  200 
kilometers  from  the  array.  The  surface  wave  train  is  complex,  with  at  least  two 
distinguishable  phases  consistently  present.  The  overall  wave  train  is  characteristic  of 
the  Lg  phase,  but  the  dominant  surface  wave  arrival  is  interpreted  as  Rg.  The  Rg  phase, 
also  referred  to  as  Rl,  and  a  later  surface  wave  arrival  (R2)  were  analyzed  to  determine 
the  effects  of  regional  and  local  geology  on  the  surface  wave  velocities  andraypaths. 
Most  of  the  data  were  analyzed  through  the  use  of  two  signal  processing  techniques, 
multiple  filter  analysis  and  cross-correlation. 
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The  surface  waves  from  mine  blasts  recorded  on  the 
Central  Minnesota  Seismic  Array  were  analyzed  for 
propagation  characteristics  and  velocities. 

The  seismic  array  was  a  six  station  vertical  array  with 
a  diameter  of  28  kilometers.  It  was  located  in  east-central 
Minnesota  and  was  emplaced  in  1976  to  monitor  local 
seismicity.  It  routinely  recorded  mine  blasts  which 
occurred  in  the  Mesabi  Range  of  northern  Minnesota,  about 
200  kilometers  from  the  array.  The  surface  wave  train  is 
complex,  with  at  least  two  distinguishable  phases 
consistently  present.  The  overall  wave  train  is 
characteristic  of  the  Lg  phase,  but  the  dominant  surface 
wave  arrival  is  interpreted  as  Rg.  The  Rg  phase,  also 
referred  to  as  Ri,  and  a  later  surface  wave  arrival  (R2) 
were  analyzed  to  determine  the  effects  of  regional  and  local 
geology  on  the  surface  wave  velocities  and  raypaths. 

Most  of  the  data  were  analyzed  through  the  use  of  two 
signal  processing  techniques,  multiple  filter  analysis  and 
cross-correlation.  Multiple  filter  analysis  was  used  to 
determine  group  arrival  time  as  a  function  of  frequency  for 
the  two  surface  waves  studied  at  each  array  station.  The 
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group  arrival  times  are  measured  at  each  frequency  after 
passing  the  signal  through  a  narrow  band  pass  filter 
centered  at  the  frequency  of  interest.  This  method  was  used 
to  calculate  group  velocities  across  the  array  by  assuming  a 
planar  wave  front  and  to  calculate  group  velocities  to  the 
Station  1  of  the  array  by  assuming  the  Rg  phase  follows  a 
straight  path  from  source  to  receiver. 

The  cross-correlation  function  is  a  measure  of  the 
similarity  of  two  signals.  It  is  sensitive  to  the  frequency 
content  and  phase  of  the  signals  being  compared  so  it  was 
used  to  identify  the  source  areas  for  the  mine  blasts,  each 
of  which  gives  rise  to  a  unique  signal  at  the  six  array 
stations.  Surface  waves  from  eight  different  source  areas 
were  distinguishable  by  cross-correlation.  Very  high 
correlation  coefficients  were  measured  at  Station  1  for 
signals  from  source  areas  separated  by  as  much  as  26 
kilometers  for  the  Rl  arrival.  This  suggests  that  the 
medium  between  the  mines  and  the  array  is  uniform.  High 
correlation  coefficients  at  Station  1  were  only  measured  for 
the  westernmost  mines. 

Phase  travel  times  were  measured  using  a  variation  of 
the  multiple  filter  method.  To  measure  phase  arrivals  at  a 
given  frequency,  the  signals  were  narrow  band  pass  filtered 
then  cross-correlated  to  determine  the  time  shift  required 
for  the  two  signals  to  be  in  phase.  Phase  differences  were 
measured  for  most  2  station  combinations  to  determine  the 
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phase  velocity  beneath  the  array  and  phase  differences  were 
measvired  at  Station  1  for  records  from  different  mines  with 
high  correlation  coefficients  to  determine  near  source  phase 
velocities. 

For  Rl,  the  group  dispersion  beneath  the  array  is 
affected  by  a  wedge-like  configuration  of  low  velocity 
Keweenawan  sedimentary  rocks  which  under ly  the  array.  The 
sequence  thickens  from  west  to  east.  High  frequencies 
propagate  at  a  nearly  constant  group  velocity  of  1.8  km/s 
across  this  structure.  At  lower  frequencies,  there  is  an 
increase  in  dispersion  observed  for  the  western  array 
stations  as  the  high  frequencies  initially  travel  at  lower 
velocities  followed  by  a  decrease  in  dispersion  as  the  low 
frequencies  are  affected  by  the  low  velocity  material.  The 
group  velocity  curves  obtained  using  all  six  array  stations 
have  a  reversed  dipsersion  trend  with  velocity  ranging  from 
approximately  1.5  km/s  at  0.4  Hz  to  1.8  km/s  at  1.5  Hz.  A 
two  layer  model  was  fit  to  the  average  group  velocity  curve 
which  requires  an  average  thickness  of  2.1  kilometers  for 
the  Keweenawan  sedimentary  rocks  beneath  the  array. 

The  phase  travel  times  for  Rl  yielded  a  nearly  constant 
phase  velocity  across  the  array  of  approximately  1.6  km/s 
for  the  range  of  frequencies  measured  (0.4-1. 2  Hz).  The 
phase  and  group  travel  times  were  also  used  to  determine  the 
azimuth  of  arrival  for  Rl.  An  average  azimuth  of  320 
degrees  was  obtained  for  the  Rl  arrival. 
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The  velocity  of  R2  across  the  array  was  about  40 
percent  lower  than  the  velocity  of  Rl.  This  suggests  that 
R2  is  a  different  mode.  The  plane  wave  approximation  did 
not  fit  the  R2  arrival  times.  The  travel  times  could  be 
matched  if  R2  is  scattered  from  a  point  source  near  the 
array.  The  scattering  has  been  interpreted  as  occurring 
near  the  contact  between  the  Keweenawan  basin  and  the 
Douglas  Fault. 

The  group  velocities  measured  to  Station  1  ranged  from 
2.65  km/s  to  2.85  km/s.  The  near  source  velocities  range 
from  1.8  km/s  to  3.6  km/s  and  confirmed  that  a  region  of 
local  high  velocity  exists  near  U.S.  Steel.  The  high 
velocity  may  indicate  the  presence  of  near  surface 
metavolcanic  rocks  within  the  Archean  basement  which 
underlies  the  Animikie  basin  near  the  source  areas. 

The  final  analysis  was  two  dimensional  raytracing  which 
attempted  to  match  the  group  travel  time  to  the  array  for  a 
1  Hz  surface  wave.  The  objective  was  to  determine  the 
velocity  gradient  required  within  the  Animikie  basin  for  the 
Rl  surface  wave  to  arrive  at  the  array  from  the  easternmost 
source  areas.  A  velocity  model  was  obtained  in  which  the 
velocity  decreases  from  2.8  km/s  to  2.0  km/s  as  the  basin 
thickens  to  the  southeast.  A  two  layer  model  with  a 
thickness  for  the  Animikie  Group  of  1.2  kilometers  was 
obtained  for  the  average  velocity  of  the  basin.  The  two 
layer  model  could  not  be  used  to  obtain  the  low  group 


vi 


velocities  implied  by  the  raytracing  model.  If  the  group 
velocities  do  include  values  as  low  as  2.0  km/s,  a  very  low 
shear  wave  velocity  must  occur  within  the  southeastern 
Animikie  basin. 
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1 . 1  Location 


The  Central  Minnesota  Seismic  Array  was  a  small 
aperture,  six  station  seismic  array  located  in  east  central 
Minnesota  (Figure  1.1).  The  geophones  were  1  Hz  analog 
vertical  instruments.  The  array  diameter  was  28  km.  Five 
stations  of  the  array  spanned  the  perimeter  of  a  circle, 
with  the  sixth  station  at  the  center  of  the  array.  The 
array  operated  from  1977  until  1982  and  was  emplaced  to 
monitor  local  seismicity  for  the  Nuclear  Regulatory 
Commission.  It  recorded  several  teleseismic  events,  a  few 
local  earthquakes,  and  routinely  recorded  quarry  blasts  at 
local  to  regional  distances  (Greenhalgh,  1979,  Mosher,  1980, 
Mooney  and  Walton,  1980) . 

Seismograms  from  blasts  from  seven  taconite  mines  of 
the  Mesabi  Range  in  northern  Minnesota  were  used  in  this 
study.  The  mines  are  large,  open  pit  quarries  which  set  off 
explosive  charges  of  several  hundred  tons.  Surface  waves 
generated  by  these  blasts  were  used  to  determine  the 
velocity  characteristics  of  the  geologic  terranes  which  lie 
between  the  source  area  (the  Mesabi  Range) ,  and  the  array 
(Figure  1.2).  A  typical  seismic  record  from  these  blasts  is 
shown  in  Figure  1.3.  There  are  several  surface  wave 


figure  l.i 


locations  in  relation  to  State  of 
Minnesota  (from  Mooney  and  Walton,  1980) 


Figure  1.2  Precambrlan  tectonic  terranes  crossed  by 

Mesabi  Range  surface  waves  recorded  on  the 
Central  Minnesota  Seismic  Array  (modified 
from  Morey  and  others,  1982,  and  Southwick 
and  others,  1988) 
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arrivals.  The  first  surface  wave  arrival,  labelled  Rl, 
closely  follows  the  S-wave  arrival  and  often  interferes  with 
the  S  coda.  Rl  is  a  Rayleigh  wave,  as  can  be  seen  by  the  90 
degree  phase  shift  between  the  vertical  and  horizontal 
components  in  the  two  component  record  of  Figure  1.4.  Rl  is 
interpreted  as  Rg.  The  dominant  frequency  of  the  surface 
waves  is  approximately  1  Hz.  There  is  other  surface  wave 
energy  present,  especially  on  Stations  2  and  5.  Some 
arrivals  have  a  higher  amplitude  on  the  horizontal 
component,  and  may  represent  Love  mode  propagation.  A 
second  surface  wave  is  consistently  present  on  several  array 
stations.  The  second  arrival,  labelled  R2  in  Figures  1.3 
and  1.4,  was  also  used  in  this  study. 

The  general  character  of  the  surface  wave  train  is  that 
of  Lg  as  described  by  several  authors.  Lg  is  a  dominant 
phase  on  regional  seismograms  (Der  and  others,  1984) .  It  is 
often  of  long  duration,  over  1  minute,  even  at  relatively 
short  distances  (Oliver  and  others,  1955) ,  as  is  the  case 
for  the  surface  wavetrain  recorded  on  the  Central  Minnesota 
Seismic  Array.  Lg  was  first  described  by  Press  and  Ewing 
(1952)  as  a  prominent  phase  on  earthquake  seismograms  with 
continental  paths.  One  of  the  most  distinctive  features  of 
Lg  is  that  as  little  as  100  kilometers  of  oceanic  crust  in 
the  propagation  path  can  extinguish  it  (Gregersen,  1984) . 

It  also  may  be  extinguished  by  large-scale  geologic  features 
(Gregersen,  1984) . 
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Lg  has  been  successfully  modelled  as  the  superposition 
of  fundamental  and  higher  Love  modes  and  higher  Rayleigh 
modes  (Kennett,  1984,  Kennett  and  Mykkelveit,  1984).  The 
modes  are  affected  by  lateral  heterogeneity  and  there  is 
much  conversion  between  modes  as  the  surface  wave 
propagates.  Kennett  (1984)  and  Kennett  and  Mykkelveit 
(1984)  modelled  and  observed  for  explosive  sources  with  an 
initial  wavefield  composed  purely  of  Rayleigh  modes,  a  net 
transfer  of  energy  from  Rayleigh  to  Love  modes  as  the  blast 
generated  surface  wave  train  propagates  through 
heterogeneous  media.  Since  the  surface  waves  observed  on 
the  CMSA  records  were  recorded  at  ne=“'  egional  distances, 
much  of  the  energy  should  '=?till  be  propagating  as  Rayleigh 
modes . 

R1  has  been  modelled  as  the  Rg  phase.  Rg  arrives  as 
part  of  the  overall  surface  wave  train  but  is  usually 
treated  separately  from  Lg  (Press  and  Ewing,  1952,  Maupin, 
1989)  .  Rg  velocity  is  controlled  by  the  velocity 
distribution  within  the  first  few  kilometers  of  the  crust 
(Maupin,  1989) .  Lg  is  sensitive  to  the  total  crustal 
velocity  distribution.  Since  the  typical  Lg  velocity  is 
3.51  km/s  and  R1  velocities  on  the  order  of  1.5  km/s  have 
been  measured,  Rl  is  most  likely  Rg. 
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1 . 2  Regional  Geology 

The  surface  waves  cross  several  major  Precambrian 
tectonic  terranes  (Figure  1.2).  The  primary  purpose  of  this 
research  was  to  describe  the  propagation  of  Rl  and  R2 
through  these  terranes.  There  are  few  exposures  of  basement 
material  in  this  region,  and  much  of  the  geology  has  been 
inferred  from  geophysical  data  (Mooney  and  others,  1970, 
Southwick  and  others,  1988,  Chandler  and  others,  1989)  and  a 
recently  completed  shallow  drilling  program  carried  out  by 
the  Minnesota  Geological  Survey  (Southwick  and  others, 

1986) . 

The  geologic  feature  which  comprises  at  least  50 
percent  of  the  total  raypath  is  the  Animikie  basin.  It  has 
recently  been  reinterpreted  as  a  foreland  basin  which  lies 
outboard  of  a  fold  and  thrust  terrane  (Figure  1.5).  The 
current  interpretation  (Southwick  and  others,  1988)  is  that 
the  basin  and  its  associated  fold  and  thrust  belt  together 
comprise  the  western  end  of  the  Penokean  orogen  and 
represent  a  series  of  early  Proterozoic  tectonic  events 
which  occurred  from  approximately  2200  Ma  to  1760  Ma.  In 
general,  the  structural  complexity  of  the  terranes 
associated  with  the  Penokean  orogen  increases  from  northwest 
to  southeast  (Morey,  1983,  Southwick  and  others,  1988).  In 
the  Mesabi  Range,  the  Animikie  Group  forms  a  southward 
dipping  sequence  of  essentially  undeformed  sedimentary  rocks 


Figure  1.5  Geologic  setting  of  the  western  Penckean 
orogen  (from  Southwick  and  others,  1988) 
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(Morey,  1983)  which,  based  on  aeromagnetic  data,  extend 
about  20  kilometers  from  the  northern  edge  of  the  basin  and 
obtain  a  maximum  thickness  of  1  km  (Chandler,  1982)  . 

The  lithologic  units  within  the  Animikie  basin  are  the 
Pokegama  Quartzite  at  the  base,  unconformably  overlying  the 
Archean  basement.  Overlying  the  Pokegama  Quartzite  is  the 
Biwabik  Iron  Formation  and  the  Virginia  Formation  (Morey, 
1983,  Southwick  and  others,  1988).  The  Virginia  Formation 
of  the  Mesabi  Range  is  correlative  with  the  Thomson 
Formation  exposed  in  the  southern  Animikie  basin.  These  are 
composed  of  intercalated  mudstone  and  siltstone  turbidite 
deposits  which  thicken  and  coarsen  progressively  from  north 
to  south  across  the  basin  (Morey  and  Ojakangas,  1970,  Morey, 
1983) .  To  the  south  of  the  Mesabi  Range,  the  structural  and 
metamorphic  grade  increases  and  as  much  as  6  kilometers  of 
sedimentary  material  may  be  present  (Ferderer,  1988) .  Much 
of  the  surface  wave  analysis  was  directed  at  determining  the 
refracting  effects  of  the  sedimentary  material  of  the 
Animikie  basin. 

A  transect  from  north-northwest  to  south-southeast 
through  the  western  Penokean  orogen  passes  through  an 
Archean  cratonic  foreland  consisting  of  a  greenstone  granite 
terrane  and  a  feature  named  the  Great  Lakes  Tectonics  Zone 
which  is  a  possible  suture  zone  between  two  Archean  terranes 
(Southwick  and  others,  1988) .  Deposited  unconformably  upon 
this  basement  material  were  sedimentary  rocks  of  the 
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Animikie  Group.  The  southern  part  of  the  basin  has  been 
folded  and  metamorphosed  so  that  tight  folds  and  a  pervasive 
regional  cleavage  were  developed  (Wright  and  others,  1970, 
Holst,  1982,  1984,  Southwick  and  others,  1988).  South  of 
the  basin  are  four  structural  panels  which  comprise  the  fold 
and  thrust  terrane  (Figure  1.5).  The  northernmost  panel,  or 
the  North  Range,  Cuyuna  District  of  Figure  1.5,  consists  of 
tightly  folded  and  thrust  faulted  pre-Animikie  volcanic  and 
sedimentary  materials.  The  southernmost  panel,  south  of  the 
Malmo  structural  discontinuity,  is  a  high  grade  metamorphic 
and  Plutonic  terrane  (Southwick  and  others,  1988)  .  Station 
1  of  the  array  was  located  over  granitic  material  of  the 
southernmost  structural  panel. 

Five  of  the  six  array  stations  were  located  over  a 
Middle  Proterozoic  basin  which  forms  part  of  the  Mid- 
Continent  Rift  System  (Figure  1.6).  The  basin  is  one  of 
sever.  I  wedge-shaped  half  grabens  which  flank  the  blocks  of 
volcanic  material  which  comprise  the  medial  portions  of  the 
rift  (Chandler  and  others,  1989) .  The  basin  probably  reaches 
a  thickness  in  excess  of  2  kilometers  in  the  array  vicinity 
(Figure  1.7)  and  the  basal  contact  may  be  faulted  in  several 
places  (Chandler  and  others,  1989) .  The  basin  fill  consists 
of  strata  of  the  Fond  du  Lac  and  Hinckley  Formations.  The 
older  Fond  du  Lac  Formation  which  comprises  most  of  the 
basin  fill  consists  of  many  fining  upward  sequences  of 
immature  elastics  (Morey  and  Ojakangas,  1983) .  The  Hinckley 
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Geologic  cross  section  along  line  B-B*  with  seismometer  locations 
superposed  (from  Mooney  and  others,  1970) 
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Sandstone  is  an  arenite  vith  98  percent  quartz  (Morey  and 
Ojakangas,  1983) . 

Chandler  and  others  (1989)  estimated  a  total  thickness 
for  the  sedimentary  section  of  the  western  basin  to  exceed  4 
kilometers  along  a  profile  (b-b*)  south  of  the  array  (Figure 
1.6).  In  addition  to  the  Hinckley  Sandstone  and  Fond  du  Lac 
Formation,  they  interpret  that  sedimentary  rocks  belonging 
to  the  Oronto  Group  are  also  present  within  the  western 
basin  along  this  profile.  Where  exposed,  the  Oronto  Group 
consists  of  conglomerate,  lithic  sandstone  and  shale 
conformably  overlying  volcanic  rocks  of  the  St.  Croix  Horst 
(Craddock,  1972).  Seismically,  the  Oronto  Group  is 
difficult  to  distinguish  from  the  Fond  du  Lac  Formation 
(Craddock,  1972,  Chandler  and  others,  1989). 

The  geometry,  structure  and  lithology  of  the  basin 
strongly  affect  the  surface  wave  propagation  across  the 
array.  The  velocity  of  the  sedimentary  material  within  the 
basin  is  much  lower  (Vp  =  3.65  km/s,  Mooney  and  others, 

1970)  than  the  velocity  of  the  other  terrenes  through  which 
the  surface  waves  propagate  (Vp  =  4. 5-6.0  km/s,  Mooney  and 
others,  1970,  Greenhalgh,  1979) .  The  low  velocity  within 
the  Keweenawan  basin  causes  the  surface  waves  to  refract  as 
they  enter  the  basin  and  the  wedge-like  geometry  affects  the 
observed  surface  wave  dispersion. 

The  basin  is  bounded  to  the  east  by  the  Douglas  Fault, 
a  steeply  dipping  northeast  trending  fault  which  separates 
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the  basin  from  the  St.  Croix  Horst  to  the  east.  The  St. 
Croix  Horst  forms  a  medial  block  of  the  rift  and  consists  of 
mainly  basalt  interbedded  with  interflow  clastic  material 
and  some  lesser  rhyolitic  material  (Morey  and  Mudrey,  1972, 
Chandler  and  others,  1989) .  There  are  also  several 
northwest  trending  faults  in  the  vicinity  of  the  array.  The 
faults  near  the  array  parallel  large  northwest  trending 
faults  which  offset  large  segments  of  the  rift  and  which  may 
be  transform  faults  (King  and  Zeitz,  1971),  or  the  northwest 
trending  faults  could  be  more  akin  to  scissor  faults 
(McSwiggen,  1987) .  Another  goal  of  this  study  was  to 
determine  which  structural  features  affect  the  surface  wave 
propagation. 

1-3  P_r_evious  Studies  for  the  Central  Minnesota  Seismic 

Arxay 

There  are  two  major  previous  studies  which  used  records 
from  the  Central  Minnesota  Seismic  Array.  Greenhalgh  (1979) 
determined  the  system  response  for  the  array  and  calibrated 
the  array  for  earthquake  detection  capability  and  array 
bias.  He  measured  the  P-wave  travel  times  to  the  array  from 
most  of  the  Mesabi  Range  mines  to  obtain  a  velocity  versus 
depth  curve  for  the  region,  and  did  some  initial  modelling 
of  the  group  dispersion  through  the  Animikie  basin.  He  used 
some  of  the  mine  blasts  to  measure  crustal  attenuation.  He 
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also  attempted  to  synthesize  the  source  waveform  for  the  P- 
waves  generated  by  these  blasts. 

Mosher  (1980)  developed  signal  processing  techniques 
for  the  array  which  were  applied  to  both  teleseismic  and 
mine  blast  data.  He  tested  several  methods  for  locating 
hypocenters.  He  described  three  Rayleigh  wave  arrivals  from 
the  mine  blasts  and  compared  several  methods  for  obtaining 
the  group  and  phase  velocities  of  the  surface  waves.  He 
used  the  velocities  he  obtained  to  interpret  the  effect  of 
the  sedimentary  wedge  beneath  the  array  on  the  surface  wave 
amplitudes,  azimuth  of  arrival,  and  dispersion 
characteristics.  He  also  did  some  preliminary  modelling  of 
the  group  travel  times  from  the  mines  to  the  array. 

1.4  Research  Objectives 

There  were  three  basic  objectives  to  this  research: 

(1)  Identify  the  source  areas  for  records  to  be  used 
in  later  analyses.  This  involved  comparing  blast 
records  released  by  the  minin*^  companies  to  our 
independently  obtained  identifications.  The 
identifications  were  verified  by  cross-correlation 
of  many  of  the  surface  wave  records. 

(2)  Identify  and  characterize  surface  wave  arrivals 
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which  can  be  correlated  on  more  than  one  station. 
This  constituted  the  majority  of  the  work  for  this 
research.  It  included  extensive  visual 
examination  of  the  seismograms  and  measurement  of 
travel  times  to  and  across  the  array. 

(3)  Describe  the  propagation  of  characteristic  surface 
waves  in  terms  of  origin,  local  and  regional 
velocity  variations  and  local  and  regional 
geology.  This  was  accomplished  by  computing  group 
and  phase  velocities  of  the  two  major  surface  wave 
'  . .vals  (R1  and  R2)  across  the  array;  measuring 
group  and  phase  velocities  for  regions  external  to 
the  array  for  the  first  surface  wave  arrival; 
computing  the  azimuth  of  arrival  for  the  two 
surface  waves;  and  two-dimensional  modelling  of 
the  raypaths  from  the  mines  to  the  array  for  the 
first  surface  wave  arrival. 

1.5  Preliminary  Work 

The  initial  stages  of  research  involved  selection  and 
digitization  of  the  original  records.  There  were  several 
hundred  mine  blast  seismograms  available.  These  were 
categorized  by  Dr.  Harold  Mooney  in  a  log  book  and  assessed 
by  him  as  to  the  quality  of  P,  S,  and  surface  waves  on  each 
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record.  From  his  log  book,  the  best  200  surface  wave 
records  were  chosen  for  digitization. 

The  records  were  digitized  on  a  large  flat  bed 
digitizer  located  at  the  Minnesota  Geological  Survey.  The 
wavefoirms  were  digitized  at  a  sampling  rate  of  3  to  4 
samples  per  half  cycle.  These  records  were  then  fitted 
using  a  cubic  spline  under  tension  (Cline,  1974) ,  and 
resampled  at  a  rate  of  0.1  seconds.  Figure  1.8  shows  an 
original  analog  record  and  the  digitized  version.  For  very 
long  duration  surface  wave  records  an  artificial  low 
frequency  component  is  apparently  introduced  by  the  cubic 
spline  routine  into  the  digital  record  (Figure  1.9).  This 
also  shows  up  quite  clearly  in  the  amplitude  spectrum  for 
this  record  (Figure  1.9),  so  the  next  step  was  to  band  pass 
filter  the  record  to  remove  the  frequency  spike  in  the 
spectrum.  The  records  were  filtered  using  a  four  point, 
zero  phase  band  pass  filter.  The  filter  was  defined  with  a 
low  cut  of  0.2  Hz  and  a  slope  of  36  db/octave.  The  high  cut 
was  4.0  Hz  with  a  slope  of  12  db/octave  (Figure  1.10). 

Figure  1.9  also  shows  the  amplitude  spectrum  after 
filtering,  and  the  resulting  time  domain  record. 

The  final  preliminary  processing  involved  designing  an 
inverse  filter  to  remove  the  instrument  response  from  the 
record.  This  is  an  important  step  for  measuring  absolute 
travel  times,  since  the  phase  response  of  the  instrument  can 
delay  a  frequency  component  by  up  to  1  second  over  the 
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Figure  1.10  Four  point  zero  phase  band  pass  filter  used  to 
remove  spline  effects 
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frequency  range  of  these  signals  (Greenhalgh,  1979) . 

Removal  of  instrument  response  is  not  a  necessary  step  when 
measuring  the  relative  arrival  time  of  a  given  frequency  at 
different  stations  since  the  instruments  presumably  have  the 
same  response  for  any  given  frequency. 

The  observed  record,  S(o)  can  be  expressed  as  a 
combination  of  the  input  into  the  system,  G((i))  and  the 
transfer  function  of  the  system,  H(c) : 

S(u)  =  G(o)  * 

To  remove  the  instrument  response,  the  spectrum  of  the 
signal  is  simply  divided  by  the  system  transfer  function. 

In  practice,  since  division  can  lead  to  instabilities,  an 
inverse  transfer  function  which  is  smoothly  varying  is 
determined,  and  the  signal  is  multiplied  by  the  inverse 
transform  function.  The  transfer  function  was  determined  by 
Mosher  (1980)  at  discrete  frequencies  and  is  listed  in  Table 
1.1.  The  inverse  phase  transfer  function  was  plotted  and 
fit  with  a  fourth  order  polynomial  (Figure  1.11).  This 
equation  was  then  used  to  remove  the  instrument  response 
after  being  tested  on  several  synthetic  records  to  assure 
that  the  correct  phase  shift  was  occurring.  An  example  of 
the  results  are  shown  in  Figure  1.12.  Overall,  the  maximum 
phase  delay  for  surface  wave  records  was  1  second. 
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Figure  1.11 


Inverse  phase  transfer  function  for  the 
Central  Minnesota  Seismic  Array 
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Table  1.1  System  Transfer  Function 


Chapter  2 

Cross-Correlation  Analysis 
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2 . 1  The  Cross-Correlation  Function 

Cross-correlation  was  chosen  as  a  method  for 
identification  of  the  blast  sources.  There  were  several 
reasons  that  the  blasts  were  unidentified.  Some  companies 
had  destroyed  records,  or  closed  down.  In  addition,  some 
mines  indicated  having  blasted  at  the  same  times.  Although 
analysis  of  relative  P-wave  arrivals  could  be  used  to 
identify  a  record,  cross-correlation  is  used  to  focus  on  the 
character  of  the  surface  wave  and  to  determine  which 
attributes  of  the  surface  waves  are  dependent  on  source 
location.  The  mining  companies  were  very  cooperative  and 
released  blast  records  for  the  period  of  operation  of  the 
array.  Many  released  maps  of  the  mines.  These  were 
especially  helpful  for  records  from  Erie  Mining  Company 
since  their  separate  pits  were  distinguishable  by  cross¬ 
correlation. 

Cross-correlation  measures  the  amount  of  similarity 
between  two  signals,  as  a  function  of  a  relative  shift  in 
time  between  the  signals.  Lee  (i960}  defines  the  cross¬ 
correlation  of  different  aperiodic  functions  as: 
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<|)l2  (t)  =  (  C)  f2(t*x)dt 

—  m 

The  correlation  function  as  shown  above  is  insensitive  to 
amplitude  differences  between  waveforms,  but  is  sensitive  to 
stretching  or  compression  of  the  waveforms  as  well  as  the 
relative  arrival  times  of  the  waveforms  being  compared. 

This,  as  applied  to  the  Mesabi  surface  waves,  makes  the 
correlation  coefficient  sensitive  to  dispersion  due  to 
differences  in  length  of  travel  path  and  differences  in  rock 
properties  along  the  travel  path.  Both  of  these  factors 
depend  on  the  source  location. 

The  most  effective  way  to  utilize  the  cross-correlation 
function  is  to  normalize  it  so  that  identical  signals  will 
have  a  maximum  cross-correlation  coefficient  of  l.  For 
digital  signals  X  and  Y  the  unnormalized  cross-correlation 
is  the  cross-product  of  X(t)  and  Y(t+T).  There  are  two 
possible  ways  to  normalize  the  cross-correlation  function 
(Neidel  and  Taner,  1971) .  The  arithmetical  normalization  is 
sensitive  to  the  amplitude  and  phase  of  the  signal.  This  is 
not  as  useful  for  correlating  signals  from  different  blasts 
in  which  effects  from  blast  size  and  ground  coupling  can 
affect  the  amplitude  of  the  signal  at  the  receiver.  What  is 
unique  for  each  blast  site  is  the  time  of  arrival  and  the 
separation  in  time  of  the  various  surface  wave  arrivals. 

The  normalization  method  called  geometrical  normalization 
(Neidel  and  Taner,  1971)  is  not  sensitive  to  amplitude,  only 
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to  frequency  and  phase.  One  result  of  using  this  type  of 
normalization  is  the  equal  weighting  of  some  long  duration, 
low  amplitude  surface  waves  with  the  primary  surface  wave 
arrivals  which  an  observer  attempting  a  visual  correlation 
would  naturally  "weight"  higher. 

These  long  duration,  low  amplitude  arrivals  correlate 
very  well  for  individual  source  areas,  and  do  not  correlate 
between  mines  at  a  given  station.  This  was  useful  in 
identifying  blasts.  The  geometrically  normalized  cross¬ 
correlation  for  two  digital  signals,  fj^  and  fj,  of  length  N 
is: 


N/2 

j‘-H/2 _ 

'  N  N 

t.  Z 

U=i  -1=1  . 

(Neidel  and  Taner,  1971) 

In  practice,  to  perform  cross-correlation  of  digital 
seismograms,  the  calculations  are  made  in  the  frequency 
domain.  Each  seismogram  station  was  compared  to  the 
corresponding  station  on  another  seismogram.  The  steps 
consisted  of: 


<^12  (i) 


1) .  Reading  both  seismograms  into  files  and  augmenting 
them  with  zeroes  so  that  each  had  the  same  2^ 
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number  of  data,  a  condition  required  by  the  Fast 
Fourier  Transform.  For  cross-correlation,  a  large 
zero  padding  is  required  to  eliminate  circularity 
(Otnes  and  Enochson,  1978) ,  so  N=2048  was  used. 

2) .  Renumbering  the  array  of  t  values  to  the  correct 
cross-correlation  offsets.  For  example,  if  one 
station  had  its  first  data  point  at  29.0  seconds 
and  the  second  at  31.0  the  cross-correlation 
function  would  have  its  first  data  point  at  t  = 
-2.0  seconds.  If  N  =  512  and  dt  =  0.1,  the 
cross-correlation  function  would  extend  from 
(-2.0  -(255*0.1))  or  -27.5  seconds  to 
(-2.0  +(256*0.1))  or  +23.6  seconds. 

Because  of  changes  in  blast  location  and 
errors  or  difficulties  identifying  the 
first  motion  on  Station  6,  the  cross¬ 
correlation  function  does  not  always 
reach  a  maximum  r-O.  For  two  sources  in 
the  same  general  location,  the  signals 
should  be  similar  enough  that  the  cross¬ 
correlation  function  peaks  very  near 
zero  offset.  In  general,  for  blasts 
from  the  same  mine,  the  cross¬ 
correlation  function  does  reach  a 
maximum  within  a  few  seconds  of  t*0.  A 
large  r  (greater  than  4  seconds)  for 
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blasts  from  the  same  mine  is  a  likely 
indication  that  different  P-wave  phases 
were  used  for  the  first  motion  on 
Station  6.  The  problem  of  using  P-waves 
as  a  reference  frame  was  a  recurring  one 
throughout  the  study.  Any  use  of  the 
total  travel  time  depends  on  the  clarity 
of  the  first  motion  of  the  seismogram. 

3) .  After  renumbering  the  t  array,  compute  the 

denominator  for  the  normalized  cross-correlation. 
This  is  constant  for  all  t,  and  can  be  expressed 
as: 

f  N  tr  U 

y  y  rj(i)2  ’ 

U=i 


Mathematically,  this  is  an  expression  for  the 
geometric  mean  of  the  energy  of  the  signals 
(Neidel  and  Taner,  1971) . 

4)  .  Reverse  time  series  1  for  cross-correlation  in  the 

frequency  domain  since  Fi*(u)F2  (<*>). 

5)  .  Compute  the  spectrum  of  each  signal  via  the  Fast 

Fourier  Transform. 

6)  .  Multiply  the  two  spectra. 
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7)  .  Take  the  inverse  transform  via  the  Fast  Fourier 

Transform. 

8)  .  Divide  each  value  of  the  function  by  the 

denominator. 

The  program  was  checked  by  computing  autocorrelations 
for  several  records.  The  autocorrelation  compares  a  signal 
with  a  time  shifted  version  of  itself.  Mathematically,  it 
is  computed  in  exactly  the  same  manner  as  cross-correlation, 
but  fi(t)  and  f2(t)  are  the  same  signal.  There  are  two 
properties  of  the  autocorrelation  function  which  make  it 
useful  for  quality  control  for  the  computer  program.  The 
autocorrelation  function,  *11(7),  is  symmetric  about  t=0 
(ie.  *11(7)  =  *11 (-7))  and  the  normalized  autocorrelation 
function  is  always  at  a  maximum  of  1  at  t=0.  Figure  2.1 
shows  an  example  of  an  autocorrelation  of  E2583,  which 
exhibits  these  properties. 


2 . 2  Cross-Correlation  for  Source  Area  Identification 


The  cross-correlation  method  worked  very  well  for 
identifying  records  from  the  same  source  area.  For  the  most 
part,  high  correlation  coefficients  (*12  >  0.7)  existed 
between  blasts  from  the  same  mines  and  low  correlations  were 
measured  for  blasts  from  different  mines.  In  general. 
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Stations  1,2,3,  and  6  consistently  have  very  good 
correlation  for  blasts  from  the  same  mines.  Station  4  and  5 
do  not  always  have  high  correlation  coefficients.  The  low 
correlation  observed  for  Station  4  can  be  explained  by  the 
fact  that  Station  4  was  very  noisy  with  many  records  only 
containing  a  few  cycles.  When  a  surface  wave  from  Station  4 
is  strong  enough  to  be  digitized,  there  is  still  a  great 
deal  of  uncertainty  about  which  surface  wave  arrival  is 
present,  and  whether  it  correlates  with  the  surface  wave 
arrivals  on  other  stations.  Since  either  different  surface 
wave  arrivals  or  only  a  few  cycles  of  a  surface  wave  are 
recorded  on  Station  4,  it  is  less  likely  that  the  records 
will  show  any  similarity. 

The  actual  identification  of  the  records  involved  more 
visual  correlation  than  actual  digital  correlation  by 
computer.  All  records  were  sorted  into  tentative 
assignments  and  all  records  from  a  single  mine  were  compared 
with  each  other  by  overlaying  them  on  a  light  table.  In 
most  cases  the  records  from  the  same  mine  would  be  in  phase 
throughout  the  record  when  overlaid.  The  match  was 
qualitatively  rated  for  each  station.  If  the  match  was 
excellent  for  all  stations,  it  was  assumed  the 
identification  was  correct.  If  the  match  was  tentative  or 
questionable  on  several  stations,  the  correlation  was 
performed  digitally.  In  some  cases,  the  digital  correlation 
was  higher  than  was  estimated  by  visual  correlation.  In 
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most  cases,  however,  the  visual  correlation  could  predict 
the  outcome  of  the  digital  correlation,  so  the  match 
between  records  was  poor,  one  of  the  records  would  be 
assigned  to  a  different  source  area.  All  ’"i.ie  blasts  were 
eventually  assigned  to  their  correct  source  area. 

If  a  blast  was  initially  misidentif ied  or  not 
identified,  it  was  correlated  visually  with  what  were 
considered  representative  blasts  from  each  source  area. 

Once  a  good  visual  match  was  obtained,  the  two  were 
correlated  digitally.  In  this  manner,  about  fifty 
unidentified  or  misidentif ied  records  were  assigned  to 
source  areas.  A  total  of  eight  source  areas  were  identified 
(Figure  1.2).  The  Butler  and  National  mines  did  not  provide 
blast  location  information,  but,  based  on  cross-correlation, 
only  one  source  location  exists  for  each  mine.  The 
distances  and  azimuths  to  these  mines  used  in  later 
calculations  are  from  Greenhalgh  (1979).  Hibbing  provided 
blast  maps  and  latitude  and  longitude  information  for  most 
of  the  Hibbing  records.  All  Hibbing  blasts  were  in  the  same 
pit,  and  have  a  very  high  maximum  cross-correlation 
coefficient,  with  little  offset  (t  =  0) ,  therefore  Hibbing 
is  considered  one  source  area  and  its  distance  and  azimuth 
are  also  taken  from  Greenhalgh  (1979) .  U.S.  Steel  blasted 
in  two  pits  during  the  time  of  the  operation  of  the  CMSA. 
These  are  called  the  east  pit  and  the  west  pit.  These 
cannot  be  distinguished  from  each  other  by  cross-correlation 
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(Figure  2.2),  therefore  U.S.  Steel  is  considered  one  source 
area  with  its  distance  to  the  array  measured  from  the  center 
of  the  west  mine  as  listed  by  Greenhalgh  (1979) . 

Erie  provided  maps  as  well  as  locations  for  each  blast. 
It  was  determined  from  correlation  results  before  seeing  the 
maps,  that  Erie  records  fell  into  two  or  more  distinct 
groups.  From  the  maps  provided  by  the  Erie  Mining  Company, 
it  was  determined  that  three  source  areas  exist  for  Erie. 
Areas  1  and  2  are  basically  different  pits,  separated  by 
about  10  kilometers.  Records  from  these  areas  do  not  show 
any  correlation  with  each  other,  but  cross-correlate  very 
well  with  records  from  their  respective  pits.  The  third 
area  for  Erie  blasts  is  to  the  northeast  of  the  main  Reserve 
pit  (Figure  2.3).  There  are  only  a  few  records  from  this 
pit,  which  is  furthest  from  the  array.  Records  from  this 
pit  do  not  correlate  consistently  with  other  records  from 
Erie  and  Reserve  (Figure  2.4).  The  Erie  pits,  as  well  as 
the  Reserve  pit,  are  in  bedrock  that  has  undergone  high 
temperature  contact  metamorphism  due  to  the  emplacement  of 
the  Duluth  Complex  (Morey,  1972).  The  metamorphism  has 
changed  the  mineral  assemblages  of  the  Biwabik  Formation  and 
other  Animikie  Group  lithologies  and  probably  causes  the  low 
correlations  observed  for  some  of  the  eastern  source  areas. 

The  Reserve  mine  consists  of  a  single  open  pit 
approximately  8  kilometers  long  (Figure  2.3).  The  company 
provided  blast  locations,  which  were  divided  into  three 
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areas.  Blast  locations  tend  to  fall  gradationally  into 
different  areas  (Table  2.1).  There  were  only  seven  or  eight 
records  from  Reserve,  so  the  records  are  treated  as  one  mine 
for  several  applications,  and  as  separate  source  areas  when 
a  change  in  distance  from  the  array  was  deemed  significant. 
The  records  from  Reserve  do  not  correlate  as  well  with  each 
other  as  records  from  other  mines.  Stations  1,  2  and  6 
generally  correlate  well.  Stations  4  and  5  (the  southern 
stations)  show  a  similarity  in  appearance  of  arrivals,  but 
there  are  small  phase  shifts  between  arrivals  on  different 
records,  i.e.,  the  relative  arrival  times  between  the  first 
and  later  wave  trains  differ  slightly.  This  can  lower 
correlation  coefficients  considerably  (Figure  2.5).  In 
addition,  blasts  from  the  southernmost  end  of  Reserve 
correlate  very  well  on  Station  1  and  2  and  occasionally  on 
Stations  3  and  6  with  records  from  Erie  area  2  (Figure  2.6). 

Table  2.1  Locations  of  Reserve  Blasts 


R826 

T59N,  R13W,  S6 

R819 

T60N,  R13W,  S32 

R2749 

T60N,  R13W,  S27 

R2157 

T60N,  R13W,  S34 

R1413 

Unknown 

R2623 

Unknown 

R1413 

Unknown 

^FTER  P  ARRIVAL  ON  STN  6.  SECONDS 
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Figure  2.6.1  Records  used  to  compare  the  correlation  between  Erie  2  and  blasts  from  the 
south  end  of  Reserve 
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This  indicates  some  comparisons  may  be  made  with  records 
from  Erie  and  Reserve. 

Other  methods  exist  for  determining  the  source  for 
these  seismic  records.  The  most  obvious,  especially  when 
using  an  array,  is  to  measure  the  relative  arrival  times  of 
P-waves  at  the  array.  These  can  be  used  to  solve  for 
azimuth  of  arrival,  assuming  a  plane  P-wave  arrives  at  the 
station.  Mosher  (1980)  deteirmined  P-wave  station 
corrections  for  the  array,  and  adapted  an  algorithm  from 
Herrmann  (1978)  to  invert  for  azimuth,  so  this  method  seems 
at  first  to  be  the  most  feasible  means  of  identifying  the 
source  areas.  This  was  not  used  initially  because  the 
cross-correlation  method  would  yield  further  information  on 
the  surface  waves.  A  P-wave  study  would  not  contribute 
significantly  to  our  knowledge  of  surface  waves.  The  cross¬ 
correlation  method  gave  an  indication  of  the  amount  of 
similarity  that  could  be  expected  of  records  from  the  same 
mine.  Also  in  the  course  of  attempting  to  identify  records, 
by  cross-correlating  records  from  different  mines,  one  could 
see  how  much  similarity  existed  between  records  from 
different  mines,  and  if  high  correlations  occur,  what 
patterns,  if  any,  can  be  found.  This  is  discussed  further 
in  the  next  section.  In  addition  Greenhalgh  (1979)  studied 
the  effect  that  errors  in  reading  P-wave  arrival  times  had 
on  the  calculated  azimuth  to  the  source.  According  to 
Greenhalgh,  if  an  arrival  time  can  only  be  read  to  the 


nearest  0.1  s,  it  will  contribute  an  error  of  1  to  2  degrees 
in  computed  azimuth.  He  also  found  the  effects  of 
wavefront  curvature  to  be  negligible  for  point  source 
distances  of  greater  than  50  kilometers,  so  a  plane  wave 
assumption  for  the  P-wave  arrivals  is  valid. 

At  a  later  point  in  the  study,  first  arrivals  were 
measured  in  order  to  estimate  origin  times  for  the  blasts. 
The  measurements  of  relative  P-wave  arrival  times  contain 
considerably  more  error  than  0.1  s  per  station.  At  these 
distances  (150  -  250  kilometers),  several  phases  of  P  may 
arrive  nearly  simultaneously  (Richter, 1958) .  The  first 
motion  at  each  station  was  measured,  and  a  quality  factor 
assigned  to  it.  No  attempts  were  made  to  match  phases, 
which  must  be  done  carefully  to  keep  from  introducing  travel 
time  error  at  these  distances.  This  was  not  an  objective  of 
the  study,  so  only  first  breaks  were  used.  These  were 
averaged  for  each  mine.  Most  stations  had  a  standard 
deviation  of  about  0.15  s.  The  azimuths  computed  from  the 
P-wave  arrivals  were  usually  correct  if  the  relative  arrival 
times  for  the  array  were  near  the  average.  A  correct 
azimuth  determination  indicates  that  identical  phases  of  P 
are  being  compared. 

The  other  method  of  identifying  a  mine,  which  was  used 
in  conjunction  with  cross-correlation,  is  the  P-wave/ surface 
wave  offset.  This  could  be  used  to  at  least  eliminate 
certain  source  areas  before  attempting  cross-correlation. 
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For  example,  if  a  surface  wave  arrived  25  seconds  after  the 
P-wave  arrived  at  Station  1  or  2 ,  the  record  was  not  from 
Erie  or  Reserve,  since  their  surface  waves  are  usually  at 
least  30  seconds  after  the  P-wave.  This  arrival  time 
corresponds  to  Butler  or  National  since  that  is  a  typical 
surface  wave  arrival  time  for  those  mines.  This  technique 
probably  could  not  be  used  to  absolutely  identify  all 
records,  since  the  onset  of  the  surface  waves  showed 
variability  within  a  source  area  and  in  some  cases  the 
difference  between  source  areas  is  only  a  few  seconds  of 
travel  time  at  most. 

2 . 3  Cross-Correlation  at  Station  1 

Although  the  primary  purpose  of  cross-correlation 
analysis  was  to  identify  source  areas,  some  geologic 
information  also  emerged.  The  most  notable  result  is  the 
high  value  of  many  correlation  coefficients  measured  at 
Station  1  for  the  Rl  arrival  from  different  mines.  The 
cross-correlation  function  at  Station  1  attains  a  value 
above  0.9  for  many  cases  where  records  from  different  mines 
are  being  compared.  This  was  first  noticed  during  attempts 
to  identify  source  areas.  The  Rl  arrival  would  correlate 
almost  perfectly  in  terms  of  alignment  of  peaks  and  troughs 
when  records  from  two  different  mines  were  overlaid  (Figure 
2.7).  In  the  case  of  Butler  and  National,  two  mines  which 
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station  1  records  from  Butler  and  National  and  their  cross-correlation.  Th 
Butler  record  has  been  shifted  4  seconds  to  show  the  alignment  of  peaks  and 
troughs . 
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are  fairly  close  together,  the  match  was  always  very  good. 
For  these  two  mines,  correlation  coefficients  as  high  as 
0.95  were  measured.  '’.e  number  of  records  with  correlation 
coefficients  above  0.9  decreases  as  the  distance  between  the 
mines  increases,  but  there  is  still  a  high  degree  of 
correlation  (*^^2  ^  0.8)  for  many  blasts  with  larger 
separations,  such  as  Butler  and  U.S.  Steel,  which  are  26 
kilometers  apart.  There  are  also  cases  where  *3^2  larger 
than  0,9  for  large  separations  between  mines.  There  are 
probably  more  for  each  mine  pair,  but  ^^.2  measured 

quantitatively  for  all  possible  pairs. 

Station  1  is  the  only  station  which  very  consistently 
has  nearly  identical  waveforms  from  different  source  areas. 
R1  at  other  stations  has  a  similar  appearance  from  different 
source  areas,  but  the  phase  match  is  generally  not  good, 
which  leads  to  lower  correlation  coefficients.  Station  1  is 
the  only  station  located  over  basement,  therefore  it  was  the 
best  station  at  which  to  study  propagation  to  the  array. 
Since  the  waves  should  not  have  propagated  through  the 
Keweenawan  sedimentary  rocks  prior  to  the  arrival  at  Station 
1,  the  measure  of  correlation  at  Station  1  can  be  used  to 
generalize  some  of  the  geologic  properties  between  the  mines 
and  the  array.  Only  the  westernmost  mines  have  high 
correlation  coefficients  measured  at  Station  1.  These  mines 
are  Butler,  Ribbing,  National,  and  U.S.  Steel.  Erie  and 
Reserve  did  not  correlate  well  with  U.S.  Steel  or  with  each 
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Other.  Even  blasts  from  different  Erie  pits  (Figure  2.8) 
did  not  have  correlation  coefficients  above  0.8. 

There  are  several  possible  ways  to  interpret  these 
results.  The  two  felt  to  be  valid  are  presented  here.  The 
most  simple  which,  without  other  measurements,  would 
probably  be  preferred,  is  to  assume  R1  takes  a  straight  path 
from  source  to  receiver.  If  this  were  the  case,  one  could 
state  that  the  western  portion  of  the  Animikie  basin  is 
homogeneous  for  the  frequency  range  monitored,  since  surface 
waves  measured  at  Station  1  from  mines  as  much  as  26 
kilometers  apart  show  no  change  of  phase  other  than  pure 
delay.  This  would  also  require  that  the  paths  through  the 
fold  and  thrust  terrane  and  the  adjacent  granitic  terrane  be 
in  homogeneous  material,  so  the  velocities  do  not  change 
with  position  within  a  terrane,  since  the  direct  path 
lengths  from  various  mines  through  these  terranes  are  fairly 
constant.  Although  these  are  not  realistic  assumptions, 
they  can  be  justified  to  some  extenc  by  the  correlation 
coefficients.  In  addition,  using  the  direct  path 
assumption,  it  can  be  stated  that  surface  waves  generated  by 
the  eastern  sources  (Erie  and  Reserve)  encounter 
heterogeneity,  as  shown  by  the  low  correlation  coefficients, 
even  for  sources  in  close  proximity.  Clearly  the  western 
Animikie  basin  is  much  more  uniform  than  the  eastern 
Animikie  basin.  The  heterogeneity  in  the  east  is  likely  due 
to  the  presence  of  Keweenawan  intr”.sx  .'“s  within  the  eastern 


Cross-correlation  of  an  Erie 
record  (E2583)  and  an  Erie  1 
record  (E84) 
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part  of  the  Animikie  basin  (Ferderer,  1988) .  Erie  and 
Reserve  are  near  the  western  boundary  of  the  Duluth  Complex, 
a  large  Keweenawan  igneous  body,  so  it  is  realistic  to 
assume  there  are  Keweenawan  igneous  rocks  within  the  eastern 
Animikie  basin,  which  surface  waves  generated  at  Erie  and 
Reserve  would  sample.  This  could  lead  to  changes  in  the 
observed  waveform  with  small  changes  in  source  location. 

A  second  explanation  for  the  high  correlation 
coefficients  measured  for  the  western  source  areas  is  that 
R1  from  all  of  these  sources  follows  a  nearly  identical  path 
from  source  to  receiver.  This  would  require  a  constant 
phase  velocity  near  the  source  since  no  additional 
dispersion  is  observed  for  the  change  in  path  length.  This 
is  more  realistic  than  assuming  the  Animikie  basin  to  be 
homogeneous  and  the  other  terranes  to  be  isotropic.  If  the 
raypaths  converge  near  the  sources,  the  dispersion  observed 
at  the  array  should  not  change  with  source  position.  This 
hypothesis  is  also  supported  by  the  measured  azimuth  of 
arrival  of  R1  at  the  array.  The  azimuth  of  arrival  at  the 
array  was  measured  by  Mosher  (1980)  using  group  arrival 
times  at  various  array  stations.  He  found  the  azimuth  of 
arrival  for  Rl  was  from  the  northwest  and  is  325  degrees  for 
all  source  areas.  The  azimuth  of  arrival  for  Rl  was 
measured  using  phase  differences  in  this  study.  The  results 
from  this  study  agree  with  Mosher's  in  that  Rl  arrives  from 
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the  northwest  for  all  source  areas.  This  is  discussed  in 
detail  ir*  Chapter  3. 

2 . 4  Conclusions 

Cross-correlation  was  a  useful  technique  for 
characterizing  the  source  area  for  both  previously 
identified  and  unidentified  records.  The  surface  wave 
arrivals  produce  unique  records  for  each  source  area. 
Cross-correlation  is  a  method  which  is  sensitive  to  the 
frequency  and  phase  of  the  records  being  compared  and  so  was 
sensitive  to  differences  in  the  records  which  occur  as  the  • 
source  location  changes.  Most  of  the  identifications  were 
made  visually,  but  were  confirmed  by  cross-correlating  the 
unidentified  record  with  an  identified  record. 

The  cross-correlation  coefficients  measured  at  Station 
1  also  yielded  information  about  the  regional  geology.  Ri 
records  from  the  western  source  areas  show  excellent 
correlation  with  each  other  on  Station  1.  This  implies  that 
the  material  being  sampled  by  Rl  prior  to  arrival  at  the 
array  is  very  uniform  with  few  changes  between  mines.  Rl 
from  the  eastern  mines  samples  more  heterogeneous  material. 


Chapter  3 


Diapersion  Analysis 


3 . 1  Introduction 

The  travel  time  as  a  function  of  frequency  is  a 
diagnostic  property  for  surface  waves.  Since  travel  time  is 
strictly  a  function  of  velocity,  if  the  source  and  receiver 
locations  are  known,  the  travel  time  as  a  function  of 
frequency  should  be  measurable.  Several  techniques  exist 
for  measuring  the  dispersion  of  a  waveform.  These  can  be 
performed  in  the  time  domain  (moving  window  analysis)  or  the 
frequency  domain  (multiple  filter  analysis,  maximum  entropy 
spectral  analysis) .  A  frequency  domain  method  called 
"Multiple  Filter  Analysis"  was  chosen  (Dziewonski  and 
others,  1969) .  It  is  efficient,  accurate,  and  still 
commonly  used  as  a  standard  method  for  processing  dispersed 
waveforms.  It  is  capable  of  separating  modes  which  overlap 
in  the  time  domain  and  can  be  used  to  calculate  group  and 
phase  velocities. 

The  multiple  filter  technique  involves  passing  the 
signal  through  a  series  of  narrow  band  pass  filters  (Figure 
3.1).  The  filters  have  different  center  frequencies  which 
correspond  to  the  frequencies  contained  in  the  signal.  For 
each  pass  band,  the  complex  envelope  is  formed  (Figure  3.2). 
Wherever  a  peak  occurs  in  the  envelope,  it  corresponds  to 
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Figure  3.1  CMSA  record.  Station  6, from  Hibbing,  top, 
and  signal  filtered  at  11  different 
frequencies 
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Figure  3. 


Envelopes  of  filtered  signals  shown  in 
Figure  3.1 
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the  arrival  time  of  a  phase  or  mode  for  the  center  frequency 
of  the  particular  pass  band.  The  results  of  the  process  are 
expressed  quantitatively  by  a  matrix  of  values  which 
represent  the*  energy  of  the  envelope  as  a  function  of  time. 
Each  column  of  the  matrix  corresponds  to  a  pass  band  center 
frequency.  Each  row  is  a  specific  time  sample  of  the 
envelope.  The  matrix  elements  are  the  energy  of  the 
envelope  normalized  to  a  maximum  of  99  db.  The  matrix  is 
contoured  and  different  surface  wave  arrivals,  and  their 
arrival  times  as  a  function  of  frequency  are  determined  by 
the  occurrence  of  peaks  in  the  envelope  matrix  (Figure  3.3). 
Each  branch  in  the  contoured  matrix  corresponds  to  a 
separate  wave  arrival.  When  applied  to  a  single  trace,  the 
multiple  filter  technique  allows  one  to  measure  the  group 
arrival  time  and  from  this,  compute  group  velocities. 

3.2  Steps  in  Multiple  Filter  Analysis 

The  flowchart  in  Figure  3.4  diagrams  the  steps  involved 
in  multiple  filter  analysis.  The  explanation  of  the  steps 
follows. 


A  seismogram  with  equally  spaced  time  samples  is 


required  for  the  use  of  a  fast  Fourier  transform.  Details 
of  the  digitization  procedure  are  covered  in  Chapter  1. 
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Figure  3.3  Amplitude  matrix  obtained  from  sampling 
signal  envelopes  at  various  center 
fre<^encies.  a  local  maximum  in  the  matrix 
indicates  the  arrival  of  a  phase  at  the 
corresponding  frequency.  Dashed  branches 
indicate  maxima  which  could  be  mistaken  as 
HI  arrivals  if  the  matrix  is  used  without 
the  seismogram. 
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Flow  chart  illustrating  steps  in  multiple  filter 
analysis 
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Step  2-Use  fast  Fourier  transform  to  get  data  into  the 
frequency  domain 

The  only  change  to  the  data  required  for  this  step  is 
to  insure  that  the  time  series  contains  2^  samples,  where  N 
is  an  integer.  This  is  accomplished  by  adding  zeroes  to  the 
trace.  Most  traces  used  a  total  of  512,  1024,  or  2048 
samples.  Only  the  longest  records,  usually  Stations  2,  5, 
and  6,  were  over  1000  time  samples  and  therefore  required 
N=2048. 

Step  3 -Selection  of  Center  Frequencies 

After  obtaining  F(q),  the  spectrum  of  f(t),  the 
function  is  multiplied  by  a  succession  of  Gaussian  filter 
functions.  The  Gaussian  filter  function  is  defined  as; 


H  ( <i))  =F  ( Q)  exp  ( -Qi 


)  ;u(3=2jrfo 


where  fg  is  the  center  frequency  of  the  pass  band  and  o  is 
the  band  width  parameter  (Dziewonski  and  others,  1969)  .  The 
range  of  center  frequencies  is  constrained  somewhat  by  the 
choice  of  a,  and  the  frequency  range  of  the  spectrum.  In 
order  to  keep  the  output  time  function  real,  all  of  the 
filter  pass  bands  must  lie  in  the  lower  half  of  the 
frequency  spectrum  which  corresponds  to  positive 


frequencies.  The  response  in  the  upper  half  of  the  spectrum 
is  constrained  by  the  condition  of  the  Fourier  transform  for 
real  time  functions  which  is;  for  N  data  points  in  the 
spectrxim,  F((i);  from  (N/2)+l  to  N  is  equal  to  the  complex 
conjugate  of  the  spectrum  from  1  to  N/2.  From  a  practical 
standpoint,  the  passband  of  the  filter  generally  extends  up 
to  about  l.Sfg.  Therefore,  for  a  frequency  sampling 
interval  df,  the  maximum  allowable  center  frequency, 

or  Fjjyg/1.3.  For  this  data,  the  time 
sampling  interval,  dt,  is  0.1  seconds.  The  Nyquist 
frequency  is  defined  as  Fj,yq=1/ (2*dt)  ,  or  5  Hz.  This  yields 
a  maximum  allowable  center  frequency  of  approximately  3.8 
Hr.  Since  the  energy  in  the  actual  signal  only  extends  to 
about  3  Hz,  this  condition  can  be  easily  met  witnout 
aliasing  problems.  In  the  actual  computations,  the  center 
frequencies  ranged  from  0.4  Hz  to  2.5  Hz  for  group  arrival 
times  and  0.3  Hz  to  1.5  Hz  for  phase  difference 
calculations. 

Step  4-Selection  of  a.  the  filter  parameter 

The  shape  of  the  filter  response  as  a  function  of  a  is 
shown  in  Figure  3.5.  The  larger  the  value  of  a,  the 
narrower  the  pass  band  and  the  higher  the  resolution  in  the 
frequency  domain.  However,  improved  resolution  in  the 
frequency  domain  causes  the  inverse  effect  in  the  time 
domain  (Dziewonski  and  others,  1969) .  In  order  to  choose  a. 


the  range  of  center  frequencies,  the  increment  of  center 
frequencies  and  the  crossover  level  of  successive  filters 
must  be  chosen  first.  For  the  group  arrival  time  analysis, 
the  center  frequencies  ranged  from  0.4  Hz  to  2.5  Hz.  1  Hz 
was  chosen  as  a  representative  frequency  from  which  to 
determine  the  filter  parameters.  A  frequency  increment  of 
0.2  Hz  gives  11  complex  envelopes  to  cover  the  frequency 
range  of  the  signal.  The  crossover  level  recommended 
between  filters  was  -12  db  (Mooney,  personal  communication) . 
The  effect  of  differing  the  crossover  level  is  shown 
schematically  in  Figure  3.6.  The  center  frequency  increment 
is  next  expressed  as  a  percentage  of  the  representative 
frequency,  which  for  this  example  is  20  percent.  In  Figure 
3.5,  the  -12  db  crossover  level  is  marked.  To  obtain  a  20 
percent  frequency  spacing,  the  filter  pass  band  and 
crossover  level  ^nould  intersect  at  0.9.  This  occurs  near 
a=50.  Since  these  parameters  were  chosen  empirically,  the 
choice  was  evaluated  by  systematically  varying  a  from  20  to 
5').  For  the  group  arrival  time  determinations,  0=50  gave 
the  best  combination  for  time  and  frequency  domain 
resolution. 

The  remaining  steps  in  the  analysis  are  automated  and 
do  not  require  any  preliminary  analysis.  Steps  5  through  8 
are  repeated  for  each  center  frequency. 


.  crossover  level  =  *25  < 


Figure  3.6  Schematic  indicating  the  relationship 

between  crossover  level  and  center  frequency 
increment .  As  the  crossover  level 
decreases,  the  center  frequency  increment 
increases .  If  the  center  frequency 
increment  becomes  too  large,  the  energy  at 
some  frequencies  will  be  lost  in  the  "aao" 
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Step  5-Window  the  Spectrum 

The  spectrum  of  the  original  signal  is  windowed  at  each 
center  frequency  by  multiplying  the  signal  by  the  Gaussian 
filter  function.  Figure  3.1  shows  an  example  of  a  filtered 
trace. 

Step  6-Obtain  the  Complex  Envelope  of  the  Windowed  Signals 
In  order  to  compute  the  instantaneous  spectral 
amplitude  of  the  windowed  spectrum,  one  must  have  knowledge 
of  the  quadrature  and  in-phase  spectrums.  The  complex 
envelope  in  the  time  domain  is  defined  as: 

C„(t)  =  An(t)  exp  (-i*n(t))  =  hn(t)  -  iqn(t) 

(Dziewonski  and  others,  1969) 

where  An(t)  is  the  instantaneous  amplitude  or  envelope  and 
q„(t)  is  minus  the  Hilbert  transform  of  h^Ct)  = 

The  definition  of  the  Hilbert  transform  (Bracewell, 1965)  is: 

Qn(«)  =  i  sgn  (u)  Hn(«) 

The  complex  envelope  can  be  determined  using  the  following 
argument,  the  quadrature  spectrum  is  defined  as: 


(1)  Qn(»)  “  i  sgn  (u)  H^  (u) 
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If  the  analytic  signal  or  complex  envelope  is; 


(2)  C^Ct)  =  h„(t)  -  iq„(t) 

and 


h(t)  f  H(w)e-^“‘=da) 

2n  J 

-CD 

00 

g(t)  =—  f  [+i  sgn  (Ci>)  e-^“*^d(o 

2n  J 

-CD 

00 

c(t)  =—  f  [1  +  sgn(o))  ]H(Q)e-^“^d(»> 

2^.4 

Equation  5  shows  that  to  obtain  the  complex  envelope  in  the 
frequency  domain,  multiply  H„((i))  by  two  for  all  positive 
frequencies  and  multiply  HnCw)  by  zero  for  all  negative 
frequencies  (Farnbach,  1975) .  The  inverse  Fourier  transform 
then  yields  the  complex  envelope  C^(t)  and  the  amplitude  may 
be  obtained  by  taking: 


(4) 

then 


An(t)  -  (Re  C„{t)2  +  im  C„(t)2)i/2 


This  is  the  envelope  as  a  function  of  time.  A  maximum  in 
the  envelope  corresponds  to  a  surface  wave  arrival  time  for 


the  center  frequency  of  the  given  pass  band.  Figure  3.2 
shows  an  envelope  obtained  from  a  filtered  trace. 


Step  7-Selection  of  the  Group  Arrival  Times 

Since  the  envelope  is  non-zero  outside  of  its  peak 
value,  the  function  is  sampled  and  contoured  in  the  time 
domain  to  locate  the  maxima  and  hence  the  group  arrival 
times.  The  envelope  may  be  sampled  at  increments  of  time 
or,  by  assuming  an  origin  time  and  distance  from  the  source, 
at  increments  of  group  velocity.  The  function  was  sampled 
at  a  constant  time  interval  of  1.5  seconds  which  allowed  the 
group  arrival  time  to  be  measured  within  0.75  seconds  out  of 
a  total  travel  time  on  the  order  of  50  to  80  seconds.  The 
accuracy  is  improved  by  averaging  the  group  travel  times  for 
each  mine.  The  arrival  times  were  later  converted  to  group 
velocities. 

Step  S-Normalization  and  Printing  of  the  Amplitude- 
Freauencv-Arrival  Time  Matrix 

Figure  3.3  shows  the  output  from  the  multiple  filter 
analysis  for  the  trace  in  Figure  3.1.  The  trace  is  also 
plotted  along  the  time  axis.  The  output  matrix  has  been 
normalized  to  99  db  for  the  peak  spectral  amplitude  for  all 
pass  bands.  This  usually  occurs  at  0.8  or  1.0  Hz.  It  is 
usually  associated  with  the  Rl  arrival,  but  it  is  not 
uncommon  for  the  peak  amplitude  to  be  associated  with  the  R2 


arrival.  The  rest  of  the  amplitudes  are  expressed  ir  db 
relative  to  the  peak  value  of  99.  If  the  actual  amplitude 
is  AMP  and  the  peak  amplitude  is  PEAK,  the  value  for  the 
pass  band  will  be: 

db  =  201ogio(AMP/PEAK)  +  99 

3 . 3  R1  Group  Dispersion  Curves  for  Array 


The  multiple  filter  technique  was  performed  on  all  of 
the  digitized  records,  using  all  channels  on  each  record. 
Many  records  were  used  to  eliminate  spurious  results  that 
could  be  present  on  any  single  channel  from  a  given  mine. 
Averaging  all  records  allowed  more  confidence  in 
interpretation  when  discussing  trends  or  calculating 
results.  Also  the  shape  of  the  dispersion  curve  can  vary 
considerably  for  a  given  source  area  and  a  given  receiver 
(Figure  3.7),  so  using  all  possible  records  allows  a 
determination  of  the  frequency  range  over  which  consistent 
results  exist.  This  is  done  by  considering  the  standard 
deviations  from  the  average  at  each  frequency. 

Figures  3.8  through  3.14  have  the  results  of  the 
multiple  filter  analysis  for  the  Rl  arrival  for  all  mines 
and  all  six  stations  of  the  array.  The  final  results  were 
obtained  after  several  preliminary  compilations.  Initially 
the  travel  times  were  measured  solely  frjm  the  amplitude 


Dispersion  curves  for  US  Steel,  Station  6.  The  curves  show  different  Rl 
dispersion  at  frequencies  above  1.2  Hz 
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grid  for  the  multiple  filter  analysis  and  the  general  shape 
of  the  dispersion  curve.  This  approach  worked  fairly  well 
for  frequencies  up  to  1.2  Hz.  Above  1.2  Hz,  there  are  many 
branches  present  in  the  seismograms  (Figures  3.3  and  3.6) 
and  it  is  difficult  to  determine  which  belongs  to  R1  without 
examining  the  seismogram.  This  preliminary  approach  led  to 
large  standard  deviations  for  the  average  arrival  time  for 
each  frequency  and  large  variations  in  group  velocity. 

The  final  compilation  involved  careful  examination  of 
the  seismograms  in  conjunction  with  the  amplitude-travel 
time  grids.  Only  those  records  with  high  quality  arrivals 
which  were  correlatable  across  the  array  were  used.  This 
procedure  often  led  to  only  two  or  three  stations  per  record 
being  used,  but  greatly  reduced  the  standard  deviations  of 
the  averages.  Using  the  seismograms  to  correlate  arrivals 
with  those  on  the  amplitude  matrix  also  increased  the 
frequency  range  over  which  the  arrivals  could  be  measured, 
since  a  branch  of  the  amplitude  grid  could  be  visually 
matched  with  the  R1  phase  on  the  record. 

A  minimum  of  five  readings  were  used  for  each  data 
point  on  the  arrival  time  curves.  If  a  station  or  frequency 
is  not  present  on  a  curve,  there  were  not  five  good  records 
available.  This  was  often  the  case  for  higher  frequencies 
and  for  Stations  4  and  5.  The  quality  of  the  data  on 
Station  4  was  generally  poor.  Station  5  had  high  quality 
data,  but  the  first  arrival  was  usually  difficult  to 
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correlate  with  the  other  stations  of  the  array,  so  often  it 
was  not  used.  The  number  of  digitized  seismograms  was  also 
a  factor.  For  example.  National  and  Hibbing  had  fifteen  or 
more  records  digitized  and  Butler  only  had  six  (since  they 
did  not  blast  often) .  There  were  many  records  available  for 
U.S.  Steel  and  Erie,  but  with  unusable  or  no  Rl  arrivals  for 
Stations  4  and  5. 

Mosher  (1980)  also  used  multiple  filter  analysis  to 
study  the  Rl  arrival.  His  observations  will  be  discussed  cn 
the  basis  of  the  group  arrival  time  curves  used  in  this 
study.  He  picked  one  record  per  mine  and  analyzed  it.  The 
results  obtained  here  differ  from  his,  because  of  the 
condition  that  Rl  be  visually  correlatable  across  the  array. 
Some  of  the  stations  he  used  had  questionable  arrivals. 

Mosher  had  the  following  conclusions  regarding  Rl 
dispersion:  Station  1  shows  the  least  dispersion,  and  the 

dispersion  increases  for  Stations  2,  5,  and  6,  then 
decreases  for  3  and  4. 

By  examining  the  dispersion  curves,  especially  the  more 
complete  ones  such  as  National  and  Hibbing  (Figures  3.9  and 
3.10),  two  things  are  apparent:  (1)  the -dispersion  is  of  the 
same  form  for  all  stations.  The  curves  have  a  lower  slope 
or  higher  amount  of  dispersion  for  the  frequency  range  0.4 
to  1.4  Hz.  Above  1.4  Hz  the  curves  dip  steeply  and  may  even 
have  a  negative  slope.  (2)  Stations  2  and  6  have  the  most 
dispersion. 
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It  is  difficult  to  discuss  the  changes  in  dispersion 
across  the  array  with  only  the  total  dispersion  curves 
available.  Since  R1  propagates  from  northwest  to  southeast 
across  the  array,  one  can  assume  that  Station  l  contains 
only  the  effects  of  propagation  through  pre-Keweenawan 
terranes.  By  removing  the  R1  dispersion  at  Station  1  from 
the  other  stations,  one  can  see  what  the  changes  in 
dispersion  are  for  each  station.  Figure  3.15  shows  the 
residual  dispersion  obtained  by  subtracting  the  group 
arrival  times  at  Station  l  from  those  at  the  other  five 
stations.  The  relative  arrival  times  are  very  consistent 
for  the  six  mines  used.  This  supports  the  assumption  that 
the  waveform  recorded  at  Station  1  is  representative  of 
propagation  to  the  array,  without  any  of  the  effects  of 
propagation  through  the  low  velocity  Keweenawan  formations. 
This  result  is  important  for  Erie  and  Reserve.  The  direct 
raypath  from  these  mines  to  Station  1  includes  this  low 
velocity  terrane,  but  the  dispersion  results  suggest  that 
the  R1  arrival  does  not  follow  a  straight  path  from  the 
mines  to  Station  l. 

A  negative  slope  on  the  composite  curves  indicates  an 
increase  in  velocity  with  frequency  and  a  positive  slope 
indicates  velocity  decreases  with  frequency.  A  vertical 
slope  indicates  no  change  in  velocity  as  a  function  of 
frequency.  A  vertical  slope  does  not  mean  that  the  velocity 
has  not  changed  between  stations,  it  can  increase  or 


FREQUENCY  (HZ)  FREQUENCY  (HZ) 


R1  COMPOSITE  CURVES  N.H.U.EI ,E2,R 


Figure  3.15  Residual  dispersion  curves  obtained  by 

removing  the  dispersion  at  Station  1  from 
the  other  array  stations 
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decrease  by  a  constant  value  for  all  frequencies. 

Stations  2  and  6  show  an  increase  in  dispersion,  or 
decreasing  velocity,  up  to  i  Hz  and  then  no  change,  or 
possibly  a  slight  increase  in  velocity  with  frequency  from 
1.2  Hz  to  2.5  Hz.  Station  4  has  reverse  dispersion  relative 
to  Station  1,  but  only  data  up  to  1.5  Hz  are  available.  The 
results  for  Station  4  are  based  on  data  from  only  National 
and  nibbing.  Stations  3  and  5  have  nearly  vertical  slopes. 
These  results  agree  with  Mosher's  with  two  exceptions.  He 
reported  an  increase  in  dispersion  for  Station  5.  He  also 
stated  that  Station  1  had  the  least  amount  of  dispersion. 
Station  1  has  approximately  the  same  dispersion  as  Stations 
3  and  5  and  Station  4  has  the  least  dispersion. 

Mosher's  interpretation  of  the  dispersion  curves  at  the 
various  stations  is  based  on  the  results  of  a  seismic 
refraction  investigation  which  characterized  the  thicknesses 
and  velocities  of  the  Keweenawan  units  in  what  would  later 
become  the  array  location.  This  investigation  found  a  low 
velocity  wedge-like  structure  beneath  the  array  (Figure 
1.7) .  The  refraction  investigation  also  found  evidence  for 
a  velocity  inversion  on  the  order  of  a  few  hundred  meters 
thickness  within  the  wedge  (Farnham,  1967,  Mooney  and 
others,  1970) .  Mosher  attributed  the  initial  increase  in 
dispersion,  which  he  observed  at  Stations  2,5,  and  6  to  the 
decreasing  velocity,  which  occurs  over  a  shorter  distance 
for  high  frequencies,  as  the  surface  waves  enter  the 
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sedimentary  wedge.  The  subsequent  decrease  in  dispersion  is 
due  to  the  velocity  inversion  which,  as  the  wedge  thickens, 
will  advance  high  frequencies  relative  to  low  frequencies. 

The  data  fit  this  interpretation,  except  at  Station  5. 
There  should  be  a  measurable  change  in  dispersion  as  the 
surface  waves  encounter  a  gradually  thickening  low  velocity 
material.  The  contact  between  the  Penokean  intrusives  and 
the  Keweenawan  sedimentary  rocks  would  have  to  be  close  to 
vertical  for  the  dispersion  to  remain  constant  at  Station  5. 

There  is  a  northwest  trending  fault  which  is  continuous 
across  the  basin.  Stations  4  and  5  are  south  of  the  fault 
(Figure  1.2).  Chandler  and  others  (1989)  describe  larger 
faults  which  this  fault  parallels  as  either  strike  slip 
faults  or  scissor  faults  which  accommodate  rift  segments  of 
alternating  symmetry.  It  is  possible  that  Stations  4  and  5 
were  located  over  a  stratigraphic  package  with  velocities 
that  are  similar  to  those  measured  by  Mooney  and  others 
(1970),  but  with  a  different  geometry  because  of  the 
presence  of  the  fault.  Qualitatively,  this  is  revealed  by 
the  magnetic  intensity  map  (Figure  3.16).  Near  Station  5, 
the  spatial  frequency  of  the  magnetic  intensity  is  high 
which  indicates  that  the  source  is  close  to  the  surface. 

The  magnetic  signal  beneath  Station  5  is  comparable  to  the 
signal  beneath  Station  1  which  is  over  basement.  This 
implies  that  the  sedimentary  section  is  quite  thin  beneath 
Station  5.  This  change  across  the  fault  could  also  account 


Figure  3.16  The  magnetic  field,  total  intensity,  near  CMSA.  Contour  interval  is 
(modified  from  Chandler,  1983  a,b) 
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for  the  difficulties  encountered  in  the  correlation  of  R1  on 
Stations  4  and  5  with  the  other  array  stations. 

3 . 4  R2  Group  Dispersion  Curves 

Group  arrival  times  for  R2  are  shown  in  Figures  3.17 
through  3.21.  The  curves  were  obtained  in  the  same  manner 
as  the  R1  dispersion  curves.  These  curves  are  more 
difficult  to  interpret  because  of  the  lack  of  uniformity 
between  sources.  However,  they  can  still  be  used  to  obtain 
group  velocities  to  compare  with  Rl.  In  general,  Station  6 
had  the  most  recognizable  R2  arrival.  It  was  usually  rhe 
highest  amplitude  arrival  on  the  record.  R2  on  Station  6  is 
characterized  by  about  three  cycles  of  signal  beyond  which 
other  arriving  energy  interferes  with  it.  Therefore  not 
much  dispersion  is  measured. 

R2  is  usually  present  on  Station  3.  The  dispersion 
curves  from  Station  3  suggest  that  high  frequencies  arrive 
earlier  than  low  frequencies.  This  may  be  an  artifact.  Low 
frequencies  may  appear  late  because  they  encounter 
interference  from  the  Rl  arrival  (Figure  3.22). 

Station  1  rarely  had  a  correlatable  R2  arrival.  No 
source  areas  had  five  useable  R2  records  for  Station  1. 

This  is  undoubtedly  due  to  R2  crossing  the  structure  of  the 
wedge  to  arrive  at  Station  1.  Station  2  also  did  not  have 
high  quality  R2  arrivals.  Only  U.S.  Steel  had  five  or  more 


group  arrival  times  for  National 


Figure  3.20  R2  group  arrival  times  for  US  Steel 


Figure  3.21  R2  group  arrival  times  for  Erie 


Figure  3.22  Group  arrival  time  amplitude  matrix  for  US 
Steel  blast  recorded  at  Station  3.  High 
amplitudes  at  frequencies  3,4,  and  5  (0.8  Hz, 
1.0  Hz,  and  1.2  Hz)  indicate  constructive 
interference  from  simultaneous  arrival  of  R1 
and  R2 . 
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records  for  Station  2.  The  dispersion  is  nearly  identical 
to  the  Station  3  dispersion  for  that  mine.  Station  4  had 
very  good  R2  arrivals,  especially  those  records  from  Hibbing 
and  National.  Early  in  the  study,  R2  on  Station  4  was  often 
mistaken  for  R1  since  it  was  the  only  surface  wave  arrival 
on  the  record.  A  very  promising  feature  on  the  curves  for 
National  and  Hibbing  is  the  similarity  of  the  R2  curve  in 
shape  and  travel  time  for  Stations  4  and  6.  If  R2  is  due  to 
reflection  or  scattering,  the  azimuth  to  the  reflector  is 
fixed  quite  well  and  does  not  change  with  frequency  for 
these  two  stations.  The  wavefront  which  defines  R2  passes 
through  these  stations  simultaneously.  It  may  be  a  plane 
wavefront  at  that  time,  or  if  a  scatterer  is  involved,  those 
stations  are  equidistant  from  the  scatterer. 

The  typical  R2  waveform  is  a  3  cycle  high  amplitude 
arrival.  It  is  very  easily  correlated  for  Stations  3,  5, 
and  6.  At  Stations  2  and  4,  it  does  not  have  the  high 
amplitudes  associated  with  3,  5,  and  6  and  appears  more 
dispersed.  It  is  rarely  correlatable  on  Station  l. 

Composite  curves  were  created  for  R2  by  taking  the 
arrival  times  of  each  frequency  relative  to  Station  3 .  The 
objective  was  not  to  obtain  changes  in  dispersion  for  R2, 
because  the  curves  varied  too  much  between  sources  for  this 
to  be  possible,  but  rather  to  look  for  consistency  in 
relative  arrival  times  at  individual  stations  for  all  source 
areas.  The  data  fall  into  fairly  broad  time  ranges  when 
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plotted  relative  to  Station  3  which  is  where  R2  arrives 
first  (Figure  3.23).  There  is  not  as  much  consistency  as 
for  the  R1  data.  Next  the  data  were  plotted  relative  to 
arrival  time  on  Station  6  (Figure  3.24).  It  can  immediately 
be  seen  that  the  curves  in  Figure  3.23  are  mirroring  the 
scatter  in  the  data  from  Station  3.  Stations  4  and  5  form 
fairly  narrow  fields  of  relative  arrival  time  when  compared 
with  Station  6.  It  appears  that  the  interference  of  Rl  and 
R2  on  Station  3  is  limiting  our  ability  to  accurately 
measure  the  arrival  time  of  either  phase  at  that  station. 
Station  3  should  probably  be  excluded  from  azimuth 
determinations  in  cases  where  it  is  not  consistent  with  the 
other  stations.  This  limits  the  possibility  of  determining 
velocity  variations  across  the  array  since  both  Stations  3 
and  4  are  unreliable.  The  data  from  Station  3  may  be  useful 
in  a  qualitative  manner  since  R2  is  consistently  present  and 
makes  a  large  contribution  to  the  signal  amplitude  at 
Station  3. 

3.5  ehgisg  Tint^.s 

The  phase  velocity  is  another  important  surface  wave 
parameter.  Two  or  more  observations  are  required  to 
determine  the  phase  velocity.  Phase  velocity  (C)  also 
varies  with  frequency  and  is  related  to  the  group  velocity 


group  arrival  times  relative  to  R2  arrival  on  Station 


R2  COMPOSITE  CURVES,  B,N,H,U,E 
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igure  3.24  R2  group  arrival  times  relative  to  R2  arrival  on  Station 
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(U)  in  the  following  manner: 

C  =*  u/k 

U  *  dw/dk  =  C  +  kdC/dk 

where  k  is  the  wavenumber  defined  as  27r/l.  Wavenumber  may 
be  dependent  on  frequency  and  position.  Phase  velocity  is 
the  actual  material  velocity,  whereas  group  velocity  is  the 
velocity  at  which  the  surface  wave  energy  propagates  (Ewing 
and  others,  1957)  so  the  azimuth  of  arrival  of  a  given 
frequency  component  is  dependent  upon  the  phase  velocity  of 
the  material.  Refraction  is  controlled  by  the  phase 
velocity,  so  that  in  order  to  measure  the  true  azimuth  of 
arrival  as  a  function  of  frequency,  knowledge  of  the  phase 
arrivals  is  required.  Also  if  accurate  phase  velocities  are 
obtained,  they  can  be  inverted  directly  for  shear  wave 
velocity  as  a  function  of  depth  (Aki  and  Richards,  1980) . 

In  quantifying  the  effects  of  the  sedimentary  wedge 
beneath  the  array  on  surface  wave  propagation,  it  was 
desirable  to  have  both  phase  and  group  velocity  information. 
The  following  quantities  were  compared  using  the  phase  and 
group  travel  time  data;  (1)  Phase  and  group  velocity 
variations  for  several  frequencies  and  array  subsets.  If 
systematic  variations  occur  in  velocity  with  respect  to 
frequency  or  position,  they  may  be  caused  by  a  change  in 
thickness  of  the  sedimentary  wedge.  The  variations  which 
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can  be  attributed  to  the  wedge  are:  (1)  a  decrease  in 
velocity  at  low  frequencies  for  the  easternmost  stations  of 
the  array,  (2)  lower  velocity  for  the  high  frequencies 
compared  to  low  frequencies  for  the  western  stations  of  the 
array  and  (3)  variation  of  measured  azimuth  for  different 
station  combinations  and  frequencies.  A  variation  of 
azimuth  should  be  measurable  for  different  array  subsets  if 
the  wedge  is  causing  continual  refraction  of  the  R1  arrival 
toward  a  direction  perpendicular  to  its  trend.  Mosher 
(1980)  reported  an  increase  in  azimuth  for  the  western 
stations,  but  he  used  incorrectly  measured  phase  lags  in  his 
phase  velocity  calculations. 

The  measurement  of  phase  velocity  is  not  straight 
forward  unless  the  station  separations  are  less  than  one 
wavelength.  For  the  array,  if  the  phase  velocity  is  between 

1  and  2  km/s  and  the  frequencies  are  approximately  0.5  Hz  to 

2  Hz,  the  wavelengths  are  on  the  order  of  0.5  kilometers  to 
4  kilometers.  The  station  separations  are  approximately  14 
kilometers,  therefore  assumptions  must  be  made  about  the 
medium  in  order  to  compute  phase  velocities. 

If  the  wavenumber  k(x,u)  does  not  vary  with  position, 
the  following  condition  holds: 

C((i>)  »  dx/dt  »  u/k(u)  (Dziewonski  and  Hales,  1972) 


If  it  is  assumed  that  all  frequency  components  have  the  same 
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initial  phase,  which  is  valid  for  explosions  (Nafe  and 
Brune,  1960) ,  the  expression  for  the  phase  velocity  becomes: 

C(u)  =  e(A  x)/A  *  =  A  x/A  t 

where  A  ♦  is  the  phase  difference  between  the  two 
seismograms  being  compared,  A  x  is  the  distance  along  the 
direction  of  propagation  between  stations  and  At  is  the 
measured  travel  time  between  stations  for  the  frequency  u. 

In  the  frequency  domain,  the  phase  difference  between 
the  seismograms  corresponds  to  the  time  domain  shift  which 
is  required  for  the  two  records  to  be  in  phase  (Bloch  and 
Hales,  1968) .  Either  of  these  may  be  determined  by 
multiplying  the  two  records.  In  the  frequency  domain,  the 
phase  of  the  cross  power  spectrum  contains  the  phase 
differences  between  the  two  seismograms  for  all  frequencies 
common  to  both  signals.  In  the  time  domain,  the  phase  shift 
between  the  two  records  can  be  determined  as  a  function  of 
time  by  cross-correlation  if  the  signal  is  a  single 
frequency. 

The  assumption  that  the  wave  number  does  not  vary  with 
position  is  a  necessary  one  f<->r  measuring  phase  velocities 
for  an  array  with  the  dimensions  of  the  CMSA.  If  the  array 
station  spacings  were  on  the  order  of  1  wavelength, 
frequency-wavenumber  analysis  could  be  used  to  determine  the 
change  in  k  with  position  and  more  accurate  phase  velocities 
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could  be  obtained.  Der  and  others  (1985)  studied  the 
spatial  coherence  and  effects  of  local  geology  on  the 
measurement  of  phase  velocities  for  the  regional  phase  Lg 
and  recommended  using  station  spacings  of  no  greater  than 
2.5  kilometers  for  determining  Lg  phase  velocities.  They 
also  found  that  the  late,  presumably  scattered  surface  wave 
arrivals  depend  strongly  on  sensor  location  and  show  almost 
no  coherence  between  stations.  Rg  and  Lg  have  similar 
frequency  content  so  the  same  conditions  should  hold  for  Rg. 
R1  and  R2  are  visually  coherent  across  the  array,  but  the 
other  surface  wave  energy  present  is  not,  so  directional 
analysis  may  only  be  performed  on  R1  and  R2 . 

In  addition,  Knopoff  and  others  (1967)  studied  the 
accuracy  of  phase  velocity  determinations  using  tripartite 
arrays.  They  found  that  the  error  in  the  measured  phase 
velocity  is  large  if  one  of  the  three  legs  of  the  array  are 
not  parallel  to  the  direction  of  the  wave  propagation.  For 
the  R1  surface  waves,  which  arrive  from  the  northwest,  this 
implies  that  the  station  pairs  1-6,  1-4,  6-4,  and  2-3  are 
the  most  reliable  for  determining  phase  velocity. 

The  method  used  to  determine  the  phase  difference 
between  two  seismograms  at  a  particular  frequency  is  a 
variation  of  the  multiple  filter  method  used  to  measure 
group  arrival  times.  The  objective  is  to  find  the  time 
shift  required  for  the  two  records  to  be  in  phase.  An  easy 
way  to  measure  this  is  cross-correlation.  The  cross- 
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correlation  function  is  the  cross  product  of  two  time  series 
measured  as  one  of  the  time  series  is  shifted  in  time 
relative  to  the  other  time  series  (see  Chapter  2  for  a 
detailed  discussion  of  cross-correlation) .  The  cross¬ 
correlation  function  is  very  sensitive  to  frequency  and 
phase.  By  applying  narrow  band  pass  filters  to  the 
seismograms  before  cross-correlating  them,  single  frequency 
records  are  created  and  the  cross-correlation  function 
becomes  dependent  only  on  phase  (Bloch  and  Hales,  1968) . 

Ideally,  the  arrival  of  interest  should  be  isolated  in 
time  from  other  modes  and  arrivals.  If  other  arrivals  are 
present,  they  are  incorporated  into  the  cross-correlation 
function  which  leads  to  multiple  branches  on  the  phase 
dispersion  curve.  For  example,  if  two  records  have  been 
filtered  about  a  center  frequency  of  1  Hz,  a  1  Hz  arrival 
corresponding  to  Rl  on  one  record  will  generate  peaks  in  the 
cross-correlation  with  the  other  record  for  both  its  Rl 
arrival  and  its  R2  arrival  (or  any  other  1  Hz  energy  present 
on  the  record) .  Therefore  it  is  necessary  to  window  the 
record  before  filtering  to  isolate  the  arrival  of  interest 
(Bloch  and  Hales,  1968). 

Rl  was  easily  windowed.  It  is  clearly  isolated  from 
other  arrivals  except  on  Station  3,  so  the  records  could 
simply  be  cut  at  some  point  between  the  end  of  the  Rl  data 
and  the  onset  of  any  later  arrivals.  The  windowing  of  R2 
was  somewhat  more  involved.  R2  is  a  large  amplitude  arrival 
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that  is  often  surrounded  by  low  amplitude  signal  (Figure 
3.25).  R2  was  isolated  in  two  steps.  First  the  record  was 
cut  on  both  sides  of  the  R2  arrival.  This  was  done  so  as  to 
include  the  complete  R2  signal  and  exclude  as  much  noise  as 
possible.  Next  a  cosine  taper  was  applied  to  both  ends  of 
the  trace.  If  a  taper  is  not  applied,  especially  if  the 
series  is  truncated  at  a  non  zero  value,  side  lobes  may 
cause  significant  high  frequency  distortion  of  the  signal 
spectrum  (Bath,  1974) .  The  tapered  window  used  was: 
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where  T  is  the  length  of  the  record,  in  seconds.  This  taper 
affects  1/10  of  the  data  on  either  enu  of  record.  The 
remaining  values  are  unchanged.  The  effect  of  the  taper  is 
also  shown  in  Figure  3.25. 

Figure  3.26  is  a  flow  chart  which  illustrates  the  steps 
involved  in  multiple  filter  analysis  for  phase  velocities. 
The  steps  are  the  same  as  those  used  in  group  velocity 
analysis,  except,  since  two  traces  are  involved,  the  cross 
power  spectrum  is  formed  (which  is  the  equivalent  of  cross¬ 
correlation  in  the  time  domain) . 
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TIME  AFTER  P  ARRIVAL  ON  STN  6.  SECONDS 


Figure  3.25  US  Steel  record,  Station  6,  before  and  after 
windowing  and  applying  a  cosine  taper  to 
isolate  the  R2  arrival 
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Figure  3.26  Flow  chart  illustrating  steps  in  multiple  filter 
analysis  for  phase  velocity  measurements 
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The  range  of  center  frequencies  used  was  0.4  Hz  to  1.8 
Hz  for  R1  and  0.3  Hz  to  1.5  Hz  for  R2.  No  attempts  were 
made  to  correlate  the  phases  above  1.8  Hz  because  there  are 
many  branches  present  on  the  amplitude  matrix  at  the  higher 
frequencies  (Figure  3.3  and  3.5).  If  more  than  one  branch 
is  present  on  the  phase  dispersion  curve,  the  mode  of 
interest  has  not  been  isolated  completely.  Since  the  range 
of  center  frequencies  was  lower  than  that  used  for  group 
velocities,  the  increment  between  center  frequencies  was 
also  lowered  to  0.15  Hz  for  Ri  and  0.10  Hz  for  R2. 

The  filter  parameter,  a,  was  again  varied  between  20 
and  50.  oe=35  had  the  best  resolution.  As  in  the  final 

group  velocity  analysis,  only  records  with  phases  which 
could  be  visually  correlated  between  stations  were  used. 

This  should  make  the  results  come  closer  to  fitting  the 
assumption  that  k(x)  is  not  varying  between  stations. 

The  cross-correlation  was  normalized  to  a  maximum  value 
of  1  as  described  in  Chapter  2.  In  addition,  in  an  attempt 
to  improve  the  results  for  R2,  the  normalization  of  the 
cross-correlation  function  was  altered  to  include 
sensitivity  to  amplitude  (Neidel  and  Taner,  1971) .  This  did 
not  lead  to  discernible  improvement  in  the  results  for  R2 
however . 

A  further  problem  arises  with  the  records  from  Station 
3.  Rl  and  R2  interfere  strongly  on  this  station  (Figure 
3.22).  Rl  masks  the  arrival  of  low  frequencies  for  R2,  and 
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the  maximum  amplitude  of  the  trace  occurs  when  Rl  and  R2 
arrive  simultaneously  and  constructively  interfere  at 
approximately  1  Hz.  Since  the  waveforms  cannot  be  isolated 
in  the  time  domain,  the  phase  differences  determined  for  R2 
are  probably  accurate  only  near  l  Hz  where  R2  has  a  large 
amplitude  pulse-like  appearance  on  the  record.  The  Rl 
measurements  are  probably  accurate  up  to  and  including  1  Hz. 

Figure  3.27  is  an  example  of  the  output  from  the 
multiple  filter  analysis  for  R2.  In  general,  the  peaks  of 
the  envelopes  were  broad,  and  the  curves  were  often 
discontinuous  between  adjacent  center  frequencies. 

Consistent  results  were  obtained  over  the  frequency  range 
0.5  Hz  to  1.3  Hz. 

The  phase  travel  time  curves  are  more  complicated  to 
display  than  the  group  travel  times  because  all  measurements 
are  made  relative  to  another  station.  For  each  frequency, 
the  number  of  measurements  made  can  be  expressed  as  a 
combination  of  six  stations  taken  two  at  a  time.  There  are 
a  total  of  fifteen  possible  phase  difference  measurements 
for  each  frequency.  In  most  cases,  fewer  than  six  stations 
were  used.  This  was  due  to  the  waveform  not  correlating  for 
certain  stations,  or  being  absent  from  a  station  altogether. 
In  some  cases,  the  phase  difference  curve  was  discontinuous 
between  adjacent  frequencies  and  was  not  used. 

Since  the  travel  times  were  used  to  compute  velocities, 
not  all  phase  differences  were  needed  except  as  quality 
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Figure  3.27  Amplitude  matrix  for  phase  velocity 

determination  for  an  R2  arrival  from  Ribbing 
recorded  on  Stations  2  and  4 
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control.  For  example,  a  three  station  velocity  calculation 
using  Stations  1,  6,  and  4  has  three  phase  differences 
available:  and  t^.g.  Since  Station  1  is  the 

first  arrival,  tg_j^  and  t^^j^  would  probably  be  used  for  the 
velocity  determination.  The  extra  measurement  t4_g  can  be 
used  to  check  the  performance  of  the  multiple  filter 
analysis  since  t^.^  should  equal  tg_j^  +  t4_g.  In  practice, 
the  travel  times  were  compiled  for  many  station  pairs  and 
those  pairs  with  the  lowest  standard  deviations  after 
averaging  were  used.  For  Rl,  the  standard  deviations  over 
most  of  the  frequency  range  were  less  than  0.5  seconds.  For 
R2,  most  of  the  standard  deviations  were  between  0.6  seconds 
and  1.2  seconds.  The  standard  deviation  for  readings  from 
Station  1  were  as  high  as  6  seconds  for  R2.  This  variation 
occurs  because  R2  is  difficult  to  identify  on  Station  1  and 
some  of  the  records  used  probably  had  a  misidentif ied  R2 
phase. 

Figures  3.28  through  3.32  are  the  phase  difference 
curves  for  Rl.  These  curves  may  be  used  to  determine  the 
changes  in  dispersion  across  the  array  and  may  be  compared 
to  the  group  dispersion  curves,  keeping  in  mind  that  the 
phase  travel  time  curves  cover  a  more  narrow  frequency  range 
than  the  group  curves. 

The  phase  travel  time  curve  for  Rl  relative  to  Station 
2  (Figure  3.29)  clearly  demonstrates  the  dispersion 
relationships  for  the  array.  The  dispersion  increases  as  Rl 


R1  PHASE  DIFFERENCES,  COMBINED  DATA 


109 


Figure  3.29  R1  phase  differences  for  Station 
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travels  from  Station  1  to  Station  2.  Stations  2  and  6  have 
approximately  the  same  dispersion.  Stations  3,  4  and  5  show 
a  decrease  in  relative  dispersion.  This  agrees  with  the 
group  dispersion  results.  One  deviation  from  the  group 
travel  time  results  is  that  the  phase  difference  curve  for 
Stations  6  and  5  (Figure  3.31)  shows  no  change  in  dispersion 
as  R1  propagates  between  these  stations.  Since  the  phase 
travel  time  curves  are  presenting  more  detail  over  a  small 
frequency  range  and  the  phase  difference  measured  does  not 
necessarily  correspond  to  the  group  travel  time,  the  trends 
can  differ.  The  phase  travel  time  curves  cover  a  range  of 
frequencies  which,  on  the  group  travel  time  curve  for 
Station  6,  shows  some  large  changes  in  slope  (Figure  3.15). 

3.6  Determination  of  Array  Group  and  Phase  Velocities 

The  travel  time  data  were  converted  to  velocities  by 
determining  azimuth  and  velocity  simultaneously.  This 
requires  travel  times  from  three  or  more  stations  in  order 
to  solve  for  two  unknowns.  The  inversion  for  azimuth  and 
velocity  assumes  a  planar  wavefront  and  a  constant  velocity 
across  the  array.  A  computer  program  from  Mosher  (1980)  for 
the  Central  Minnesota  Seismic  Array  was  used.  The  program 
performs  a  least  squares  inversion  of  the  travel  times  to 
solve  for  azimuth  and  apparent  velocity.  The  program  also 
computes  standard  errors  in  azimuth  and  velocity  and  station 
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residuals.  The  station  residuals  were  useful  because  they 
reveal  which  travel  times  are  most  difficult  to  fit  with  the 
plane  wave,  constant  velocity  model. 

The  group  and  phase  velocities  of  Rl  and  R2  were 
calculated  using  five  or  six  stations  for  all  of  the  mines. 
In  addition,  the  velocities  were  determined  for  all  possible 
three  station  combinations.  This  was  done  to  look  for 
variations  in  azimuth  and  velocity  which  are  station 
dependent  and  attributable  to  the  wedge  structure. 

The  anticipated  results,  based  on  the  travel  time 
curves  were: 

(1)  A  decrease  in  azimuth  with  increasing  frequency 
for  Rl.  Upon  entering  the  wedge,  the  high 
frequencies  should  experience  a  higher  initial 
velocity  gradient  and  be  refracted  more  strongly. 
This  could  also  manifest  as  a  higher  azimuth, 
especially  at  low  frequency,  for  the  westernmost 
stations  of  the  array  (1,2, 5, 6). 

The  western  stations  of  the  array  should  also  have 
higher  velocities  at  low  frequencies  when  compared 
to  the  eastern  stations.  The  possible  presence  of 
a  velocity  inversion  makes  it  difficult  to  predict 
velocity  trends  for  higher  frequencies. 


(2) 
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3.6.1  Previous  Results 

Mosher  (1980)  reported  a  decrease  in  azimuth  with 
increasing  frequency  for  a  single  record  from  Reserve.  He 
also  reported  phase  velocities  for  R1  which  vary  from  1.5 
km/s  to  2.0  km/s  as  an  average  for  the  whole  array  using 
individual  blasts  from  each  mine.  Mosher  defined  a  shallow 
subset  of  the  array  as  Stations  1,2, 5, 6  which  had  a  phase 
velocity  range  from  1.7  km/s  to  2.0  km/s  and  a  deep  subset 
as  Stations  2, 3, 4, 6  which  had  a  constant  phase  velocity  of 
1.4  km/s. 

Figure  3.33  shows  a  simple  two  dimensional  wedge  and 
the  group  velocity  curves  associated  with  the  structure. 
These  were  used  by  Mosher  to  qualitatively  explain  the 
changes  in  dispersion  which  occur  as  the  wedge  thickens. 

The  form  of  the  curves  roughly  resembles  the  group  travel 
time  curves,  especially  for  Stations  2  and  6.  There  is  a 
pronounced  group  velocity  minimum  for  1  Hz  for  a  wedge 
thickness  of  about  0.9  km.  Although  no  array  subgroups  had 
R1  velocities  as  low  as  1.2  km/s  for  0.8-1  Hz  (the  dominant 
frequency  range  of  the  signal) ,  Station  2  has  high  amplitude 
R1  arrivals  and  Station  6  has  the  highest  R2  amplitudes  so 
these  stations  may  be  located  near  a  group  velocity  minimum 
for  the  respective  arrivals  (Ewing  and  others,  1957) . 
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3 . 7  Velocity  Results  for  R1 

The  group  velocities  for  the  array  as  determined  from 
R1  travel  times  are  summarized  in  Table  3.1  and  shown  in 
Figure  3.34.  The  general  trend  is  toward  increasing 
velocity  with  increasing  frequency.  The  plane  wave 
assumption  worked  well,  especially  for  frequencies  near  1 
Hz.  The  results  for  U.S.  Steel  were  somewhat  poor,  but 
included  only  four  stations. 

There  were  not  easily  discernible  trends  in  group 
velocity  or  azimuth  for  array  subsets.  The  westernmost 
station  groups  1,2,5  and  1,5,6  have  the  predicted  higher 
group  velocities  and  somewhat  higher  azimuths,  especially 
for  the  lower  frequencies  such  as  0.4  Hz  and  0.8  Hz  (Table 
3.2).  The  increase  in  velocity  is  fairly  large  compared  to 
the  six  station  results,  about  0.25-0.50  km/s.  The 
azimuths  increase  about  10  degrees  over  all  frequencies. 

For  comparison  purposes,  group  velocity  results  for  two 
frequencies,  0.8  Hz  and  1.2  Hz,  are  presented  along  with 
interpretation  of  the  results. 

The  set  of  group  velocities  for  0.8  Hz  consisted  of  45 
three  station  velocity  and  azimuth  determinations.  All 
source  areas  except  Reserve  were  used.  The  velocities  fall 
into  two  groups;  velocities  greater  than  1.5  km/s  and 
velocities  less  than  1.5  km/s.  For  Rl,  normal  azimuths  were 
defined  as  any  value  between  305  degrees  and  330  degrees. 
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Figure  3.34  R1  group  velocity  as  a  function  of  frequency,  determined  using  four  or  more 
stations  for  each  source  area 


Table  3.1  Rl  Group  Velocity  and  Azimuth  Results  Across  CMSA 

(Velocities  are  given  in  km/s,  azimuths  in  degrees) 
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Table  3.2  R1  Azimuth  and  Group  Velocity  for 
Station  1/2/5  and  1/5/6 
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Of  33  group  velocities  over  1.5  km/s,  28  had  normal 
azimuths.  These  were  averaged  to  give  an  array  group 
velocity  of  1.73  km/s  and  a  group  azimuth  of  319  degrees. 
Azimuths  which  were  high  usually  had  Stations  1  or  2  as  one 
of  the  three  stations.  These  azimuths  range  from  323 
degrees  to  332  degrees  and  indicate  a  change  in  the  azimuth 
of  the  wave  packet  as  it  crosses  the  array  from  west  to 
east.  These  higher  azimuths  occurred  both  for  fast  station 
groups  (1,2,5;  1,2,4;  1,2,6)  and  slower  station  groups 
(1,2,3) . 

There  were  12  group  velocities  for  0.8  Hz  which  were 
less  than  1.5  km/s.  The  average  is  1.36  km/s,  0.4  km/s 
slower  than  the  main  group.  The  azimuth  is  the  same.  These 
three  station  groupings  all  have  Station  2  in  common,  but 
not  all  data  from  Station  2  fall  into  this  low  velocity 
group.  The  case  for  a  group  velocity  low  near  Station  2  is 
strengthened  by  the  high  amplitudes  observed  at  this  station 
(Ewing  and  others,  1957).  Mosher's  theoretical  group 
velocity  curves  (Figure  3.33)  for  a  wedge  structure  indicate 
low  group  velocities  occur  at  several  frequencies  as  the 
wedge  thickens  between  approximately  0. 6-1.0  kilometers. 

The  group  velocities  (U)  for  1.2  Hz  fall  into  three 
categories;  U  <  1.5  km/s,  1.5  <  U  <  2.0  km/s,  and  U  >  2.0 
km/s.  There  were  a  total  of  39  velocity  determinations  for 
this  frequency.  Of  these,  7  are  less  than  1.5  km/s,  23  are 
between  1.5  km/s  and  2.0  km/s,  and  9  are  greater  than  2.0 
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km/s.  For  the  low  velocity  group,  6  of  the  7  have  Station  2 
as  one  of  the  three  stations. 

Station  4  is  associated  with  the  high  velocity  group 
for  7  of  the  9  high  velocities,  therefore,  there  may  be  a 
velocity  high  associated  with  Station  4  at  higher 
frequencies.  Since  the  data  quality  was  generally  poor  for 
this  station,  and  the  reliability  of  the  travel  time 
readings  decreased  as  frequency  increased,  this  hypothesis 
needs  validation  from  either  the  phase  velocity  results  or 
the  R2  results. 

The  remaining  group  velocities  for  1.2  Hz  lie  between 
1.5  km/s  and  2.0  km/s.  Of  these  23,  18  have  normal  azimuths. 
The  average  group  velocity  is  1.71  +  0.24  km/s,  nearly 
identical  to  the  group  velocity  for  0.8  Hz.  The  average 
azimuth  is  322  ±  7  degrees,  essentially  the  same  as  the 
result  for  0.8  Hz. 

A  theoretical  group  velocity  versus  period  curve  has 
been  fit  to  the  average  array  group  velocities  (Figure 
3.35).  The  curve  assumes  a  two  layer  model  and  was  computed 
using  a  program  from  Wang  (1985).  The  velocities  and 
densities  for  the  model  are  shown  in  Figure  3.36.  The 
velocities  were  taken  from  Mooney  and  others  (1970) .  The 
densities  were  taken  from  Chandler  and  others  (1989).  The 
best  fit  to  the  data  was  obtained  using  a  thickness  for  the 
Keweenawan  sedimentary  layer  of  2.1  kilometers.  This  could 
represent  the  average  velocity  structure  beneath  the  array. 
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Figure  3.36  2  layer  model  used  to  fit  the  group 
velocities  measured  for  CMSA 
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The  fit  to  the  observed  values  is  excellent.  The  dispersion 
curves  for  R1  (Figure  3.15)  presented  earlier  in  this 
chapter  can  also  be  interpreted  using  this  model  by 
comparing  the  changes  which  occur  for  the  theoretical  group 
velocity  as  the  sedimentary  layer  thickens.  At  very  high 
frequencies,  the  group  velocity  will  remain  constant.  This 
explains  why  the  dispersion  curves  are  vertical  at 
frequencies  higher  than  approximately  1.5  Hz.  As  the  wedge 
thickens,  a  group  velocity  minimum  occurs  which  migrates 
toward  lower  frequencies.  Initially,  the  dispersion  will 
increase  since  low  frequencies  will  travel  much  faster  than 
high  frequencies.  As  the  group  velocity  minimum  passes 
progressively  through  0.8  Hz,  0.6  Hz  and  0.4  Hz,  the 
dispersion  will  reverse  as  the  high  frequencies  propagate  at 
a  constant  and  higher  velocity  than  the  lower  frequencies. 
This  model  does  not  require  a  velocity  inversion  to  explain 
the  dispersion  changes  across  the  array.  The  model  does 
imply  that  initially  the  wedge  thickens  quite  rapidly  since 
no  velocities  as  high  as  those  for  the  h=0.75  kilometers 
curve  are  observed. 

An  identical  analysis  was  performed  on  the  phase  travel 
times  for  R1 .  The  phase  velocity  as  a  function  of  frequency 
is  presented  in  Figure  3.37  and  Table  3.3.  The  plane  wave 
approximation  worked  very  well.  The  standard  errors  in 
azimuth  were  less  than  5  degrees  for  most  frequencies  and 
the  errors  in  velocity  were  on  the  order  of  .05  km/s.  The 


126 


(S/WX)  AiI0013A  3SVHd 


Figure  3.37  R1  phase  velocities  obtained  for  the  CMSA.  N,H,  and  E  had  five  or  more 
phase  difference  determinations  per  frequency.  The  combined  data  curve 
includes  all  eight  source  areas 
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phase  velocities  range  from  1.63  km/s  to  1.72  km/s,  which  is 
less  range  than  the  group  velocities  have. 


Table  3.3  Rl  Phase  Velocity  and  Azimuth/ 
Results  for  All  Source  Areas  Combined 


Frequency  (Hz) _ Phase  Velocity  (km/s)  Azimuth  (Degrees) 


l.ss"""" . . 

1.626  ±  .102 

318  ±  6 

CO 

• 

1.628  ±  .087 

319  ±  5 

CM 

00 

• 

1.664  ±  .037 

321  ±  2 

.96 

1.762  ±  .052 

319  ±  2 

1.09 

1.795  ±  .051 

319  ±  3 

1.23 

1.716  ±  .063 

320  ±  4 

The  phase  velocities  were  also  computed  for  most  three 
station  combinations.  The  azimuths  and  velocities  were  very 
consistent  and  have  a  slight  trend  toward  increasing 
velocity  with  increasing  frequency.  The  average  velocity 
for  the  three  station  combinations  increases  from  1.51  km/s 
at  0.55  Hz  to  1.69  km/s  at  1.23  Hz. 

Station  4  tends  to  have  high  phase  velocities 
associated  with  it  at  frequencies  of  1.0  Hz,  l.i  Hz  and  1.2 
Hz.  The  phase  velocities  determined  using  Station  4  for 
these  frequencies  were  close  to  2  km/s.  This  result 
confirms  the  group  velocity  results  which  suggest  a  high 
frequency  velocity  high  in  the  vicinity  of  Station  4.  A 
high  frequency  anomaly  implies  an  increase  in  the  velocity 
of  the  sedimentary  layer,  possibly  due  to  the  presence  of 
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volcanic  material  in  the  detritus  as  the  rift  is  approached. 

The  station  combination  1,2,5  also  has  higher  than 
average  phase  velocities.  The  phase  velocity  for  this 
station  combination  is  frequency  dependent.  The  average 
value  over  all  frequencies  is  1.85  km/s.  The  velocity  is 
2.0  km/s  for  0.55  Hz  and  0.7  Hz,  but  drops  to  about  1.7 
km/s,  the  average  for  the  rest  of  the  array,  for  frequencies 
higher  than  0.7  Hz.  Phase  velocities  were  not  determined 
for  the  combination  1,5,6. 

Since  no  changes  in  azimuth  as  a  function  of  stations 
used  were  detected  for  R1  using  phase  travel  times,  it 
appears  that  R1  actually  does  arrive  from  the  northwest.  A 
second  possibility  is  that  since  the  velocities  beneath  the 
array  are  so  low,  most  of  the  refraction  of  the  Rl  arrival 
as  it  enters  the  sedimentary  wedge  occurs  over  a  short 
distance  and  cannot  be  detected  by  the  array.  This  is 
addressed  further  in  Chapters  4  and  5. 

3.8  Summary  of  Array  Azimuth  and  Velocity  for  Rl 

The  velocity  and  azimuth  results  for  Rl  may  be 
summarized  as  follows: 

(1)  Both  group  and  phase  arrivals  for  any  frequency 
may  be  approximated  by  a  plane  wave  which 
propagates  across  the  array  from  northwest  to 
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southeast  at  a  nearly  constant  velocity  of  1.7 
km/s  with  an  azimuth  of  320  degrees. 

(2)  There  is  a  group  velocity  low  associated  with 
Station  2 .  This  was  also  described  by  Mosher 
(1980),  and  should  explain  the  high  amplitudes 
associated  with  the  R1  arrival  at  Station  2 . 

(3)  The  average  group  velocity  structure  beneath  the 
array  can  be  modelled  as  a  two  layer  structure 
which  is  interpreted  as  Keweenawan  sedimentary 
rocks  over  basement.  The  required  thickness  of 
the  sedimentary  material  is  2.1  kilometers. 

(4)  There  are  higher  group  and  phase  velocities 
associated  with  the  station  combinations  1,2,5  and 
1,5,6.  The  phase  velocity  high  occurs  only  for 
low  frequencies.  The  high  velocities  are 
associated  with  the  sedimentary  materials  being 
thin  beneath  Station  5  and  nonexistent  beneath 
Station  1.  There  is  also  a  high  velocity 
associated  with  Station  4  at  frequencies  of  i  Hz 
or  higher. 

(5)  There  is  group  velocity  evidence  of  changes  in 
azimuth  as  R1  propagates  across  the  array.  There 


is  no  evidence  from  the  phase  velocities.  The 
errors  in  azimuth  are  small  when  all  six  stations 
are  inverted  for  azimuth  and  phase  velocity  for 
both  data  sets.  This  implies  that  the  curvature 
of  the  wavefront  is  small,  and  the  orientation  of 
the  wavefront  does  not  change  as  Rl  propagates 
across  the  array.  The  velocity  beneath  the  array 
is  so  low,  that  any  wavefront  with  non-normal 
incidence  will  undergo  considerable  refraction  as 
it  enters  the  basin.  Changes  in  the  angle  of 
refraction  may  not  be  spatially  resolvable  by  this 


array. 


The  velocity  analysis  for  R2  proceeded  in  the  same 
manner  as  for  Rl.  The  travel  times  for  the  six  stations 
were  inverted  assuming  R2  is  planar.  It  was  anticipated 
that  R2  would  have  velocities  close  to  the  Rl  velocities. 
This  was  not  the  case.  R2  has  a  lower  velocity  than  Rl. 
Another  complication  for  R2  is  that  its  azimuth  and  velocity 
are  strongly  dependent  on  the  station  combination  used  to 
determine  them.  This  implies  that  R2  is  not  planar. 

The  velocities  for  R2  are  presented  in  Figure  3.38  and 
Table  3.4.  The  group  velocities  for  National  and  U.S.  Steel 
are  based  on  six  stations.  The  curve  from  Butler  is  based 
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Figure  3.38  R2  group  velocity  as  a  function  of  frequency  obtained  using  four  or  mor 
array  stations  per  mine.  R2  phase  velocity  as  a  function  of  frequency 
obtained  using  data  from  all  source  areas. 
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on  five  stations  and  does  not  include  data  from  Station  1. 
The  velocities  are  in  good  agreement  with  each  other.  The 
velocities  range  from  1.05  km/s  to  1.60  km/s,  considerably 
lower  than  the  group  velocity  of  Rl.  The  standard  errors  in 
velocity  are  fairly  low  (.04  km/s  to  .13  km/s). 

The  azimuths  are  internally  consistent.  National  and 
U.S.  Steel  are  also  consistent  with  each  other.  The  azimuth 
for  Butler  is  about  10  degrees  higher  than  the  other 
azimuths.  This  is  related  to  the  absence  of  Station  1  and 
is  discussed  further  in  the  next  section. 


Table  3.4  R2  Group  Velocity  and  Azimuth  Results 


Frequency  (Hz)  B~Veloeity _ N»Velocity _ D»Velocity 


0.4 

.  — 

1.15  ±  .09 

0.6 

1.36  ±  .35 

1.20  ±  .13 

1.30  ±  .06 

to 

• 

o 

1.17  ±  .11 

1.16  ±  .08 

1.26  ±  .06 

1.0 

1.23  ±  .10 

1.14  ±  .05 

1.23  ±  .05 

1.6 

1.10  ±  .04 

1.105  ±  .06 

1.24  ±  .06 

Frequency  (Hz) _ B-Azimuth _ N~Azimuth _ U~Azimuth 


0.4 

68  ±  7 

0.6 

77  ±  8 

68  ±  10 

68  ±  4 

CO 

• 

o 

78  ±  9 

63  ±  6 

64  ±  4 

1.0 

75  ±  9 

64  ±  4 

66  ±  4 

1.6 

73  ±  3 

66  ±  5 
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The  phase  velocities  for  R2,  which  were  determined  from 
data  combined  from  several  source  areas,  are  shown  along 
with  the  group  velocities  from  R2  in  Figure  3.38.  The  phase 
velocities  are  similar  to  the  group  velocities.  The  errors 
in  phase  velocity  are  small,  less  than  0.07  km/s.  The 
azimuths  are  also  similar  to  those  obtained  from  group 
travel  times.  The  six  station  phase  velocity  results  for  R2 
appear  to  support  the  assumption  that  R2  may  be  approximated 
by  a  plane  wave  which  propagates  across  the  array  at  a 
constant  velocity.  The  three  station  velocities  yielded 
much  different  results. 

3.9.1  Velocity  Results  for  Three.  Four  and  Five  Stations 

The  more  easterly  azimuths  which  were  determined  from 
the  R2  group  times  from  Butler  suggest  that  Station  l  has  a 
significant  effect  on  the  computed  azimuths  and  velocities 
of  R2.  To  determine  how  the  azimuth  and  velocity  are 
affected  by  a  single  station,  three,  four  and  five  station 
azimuths  and  velocities  were  computed  for  R2.  The  azimuth 
of  R2  varies  considerably  depending  on  which  stations  are 
used. 

The  two  stations  which  most  strongly  affect  the  azimuth 
when  using  four  or  more  stations  are  1  and  2.  The  other 
four  stations,  when  taken  in  combinations  of  three  and  four, 
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have  azimuths  and  velocities  nearer  to  the  six  station 
averages . 

The  group  travel  times  from  National  were  studied, 
since  all  six  stations  were  available.  When  Station  2  is 
eliminated  the  azimuth  decreases  (at  0.8  Hz  for  example,  the 
azimuth  decreases  from  63  degrees  to  55  degrees) .  When 
Station  2  is  included  and  Station  1  is  excluded,  the  azimuth 
increases  to  68  degrees  (for  0.8  Hz).  There  is  no  change  in 
velocity  in  either  case.  This  trend  holds  for  all 
frequencies  and  also  agrees  with  the  azimuth  results  from 
Butler.  When  Stations  1  and  2  are  eliminated,  the  azimuth 
and  velocity,  for  all  frequencies,  is  about  60  degrees  and 
1.25  km/s.  Since,  when  Station  2  is  included,  the  azimuth 
increases  to  a  more  easterly  direction,  R2  at  Station  2  is 
late  relative  to  Station  3,  the  first  arrival.  The  solution 
increases  the  distance  the  wavefront  travels  to  Station  2 
relative  to  Station  3  (Figure  3.39).  Likewise,  the  low 
azimuth  associated  with  Station  1  implies  that  it  is 
relatively  early.  This  supports  the  general  trends  observed 
for  Rl.  However,  there  are  more  complications  for  R2. 

When  the  travel  times  for  R2  were  used  in  three  station 
combinations,  there  were  several  anomalous  results.  Normal 
azimuths  were  defined  as  any  value  between  50  and  75 
degrees.  The  range  of  normal  velocities  was  defined  as  i.O 
km/s  to  1.5  km/s.  There  were  four  combinations  that  had 


Figure  3.39  Schematic  illustrating  how  the  plane  wave 
solution  for  R2  travel  times  will  obtain  a 
more  easterly  azimuth  if  Station  1  is 
excluded  and  Station  2  is  relatively  late. 

The  lines  to  the  wavefront  represent  relative 
travel  time  from  Station  3.  Station  2  has 
the  largest  change  in  relative  arrival  time 
for  the  two  azimuths  shown. 
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systematic  deviatiations  from  the  average  results.  These 
are  listed  in  Table  3.5. 

The  ranges  of  azimuth  and  velocity  for  R2  suggest  that 
it  is  not  planar.  It  is  a  coherent  arrival,  especially  on 
Stations  3,  5,  and  6.  To  explain  the  deviations  in  azimuth 
and  velocity,  it  was  modelled  as  surface  wave  emanating  from 
a  point  source,  such  as  a  diffractor,  located  near  the  array 
in  the  vicinity  of  the  Douglas  fault  (Figure  3.40).  This 
type  of  source  can  appear  planar,  especially  when 
considering  only  the  southern  stations  of  the  array  (Figure 
3.41)  . 

The  velocity  from  the  point  source  to  the  receivers  was 
assigned  a  constant  value  of  1.2  km/s.  The  computed  travel 
times,  and  actual  travel  times  relative  to  Station  3  are 
included  in  Figure  3.40.  The  location  of  the  source  matches 
the  phase  travel  times  for  0.9  Hz.  These  were  the  most 
consistent  travel  times  measured  for  R2. 

The  azimuth  and  velocities  obtained  for  the  model  are 
listed  in  Table  3.5.  The  azimuth  obtained  to  the  point 
source  and  the  velocity  across  the  array  for  the  six  station 
inversion  essentially  match  the  results  obtained  when  a 
plane  wavefront  is  assumed.  The  azimuths  and  velocities  of 
the  three  station  groups  show  the  same  trends  exhibited  by 


the  real  data. 


Tabla  3.5  Anomalous  Azimuth  and  Velocity  Results  for  R2 

A.  Measured  from  Travel  Time  Data  (All  frequencies) 
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Stations 

Azimuth  (Degrees) 

Velocity  (Icm/s) 

1,2, 3, 4, 5. 6 

69  ±  5 

1.20  ±  .12 

3,5,6 

49 

1.2 

55 

2.0 

2,5,6 

80 

1.0 

1.3,6 

92 

1.3 

2,3,6 

76 

1.5 

B.  Computed  from  Point  Source  Model  for  0.9  Hz 


1,2, 3, 4, 5, 6 

60  ±  6 

ibkhihi 

3.5,6 

47 

1.19 

wBommam 

27 

2.40 

2,5,6 

78 

1.04 

1,3,6 

87 

1.48 

2,3,6 

68 

1.43 

C.  Computed  from  Point  Source  Model  with  Low  Velocity  at 
Station  2 


1,2, 3, 4, 5, 6 

63  ±  6 

1.20  ±  .08 

2,4,6 

44 

2.85 

2,5,6 

86 

.96 

2,3,6 

76 

1.47 
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Since  there  is  evidence  for  a  low  velocity  associated 
with  Station  2,  the  model  was  modified.  The  velocity  from 
the  source  to  Station  2  was  decreased  from  1.2  km/s  to  l.i 
km/s.  The  results  are  also  in  Table  3.5.  The  agreement 
between  this  model  and  the  actual  data  is  excellent,  except 
for  the  station  combination  2,4,6.  The  travel  time  to 
Station  4  is  off  by  2.8  seconds.  There  is  evidence  from  Rl 
of  a  velocity  high  associated  with  Station  4  at  this 
frequency.  It  is  not  necessary  to  invoke  an  increase  in 
velocity  near  Station  4  because  the  station  combination 
2,4,6  is  not  geometrically  reliable  for  R2.  The  procedure 
was  attempting  to  fit  small  differences  in  travel  times  to  a 
widely  spaced  station  group  (Figure  3.42).  Velocities 
measured  perpendicular  to  the  wavefront  are  considered  more 
reliable  than  measurements  made  parallel  to  the  wavefront 
(Knopoff  and  others,  1967) . 

The  last  consideration  for  R2  is  its  low  velocity  when 
compared  to  Rl.  Kennett  (1984)  found  that  for  explosive 
sources,  which  have  an  initial  wave  field  of  purely  Rayleigh 
modes,  as  the  surface  wavetrain  progresses  through  a 
heterogeneous  medium,  there  is  a  net  tremsfer  of  modes  from 
Rayleigh  to  Love  modes.  This  cannot  be  verified  without 
three  component  data  for  which  the  radial  and  tranverse 
components  may  be  isolated,  but  is  consistent  with  the 
complete  surface  wave  train  being  Lg  which  is  a  combination 
of  many  Rayleigh  and  Love  modes.  R2  arrives  from  the 
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10  km 


Figure  3.42  Schniematic  illustrating  how  a  large  change  in 

wavefront  orientation  does  not  result  in  a  large 
change  in  relative  arrival  times  for  station 
combinations  which  are  nearly  parallel  to  the 
wavefront . 
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northeast  and  has  a  large  east-west  horizontal  component. 
This  implies  that  there  is  a  significant  SH,  or  Love  mode 
contribution  to  the  amplitude  of  R2.  If  R2  is  a  converted 
mode,  it  could"  travel  at  a  different  velocity  from  Ri. 

3.9.2  Summary  of  Azimuth  and  Velocity  Results  for  R2 

The  velocity  and  azimuth  results  for  R2  may  be 
summarized  as  follows: 

(1)  R2  is  generated  near  the  Douglas  fault  and  may  be 
considered  a  scattered  wave.  The  azimuth  to  the 
point  source  from  the  center  of  the 

array  is  approximately  65  degrees. 

(2)  The  velocity  of  R2  is  1.2  km/s  and  is 
constant  over  all  frequencies.  There  is  a 
low  velocity  associated  with  R2  at  Station  2 . 

(3)  R2  may  not  be  the  same  surface  wave  mode  as 
Rl,  therefore  it  cannot  be  used  to  determine 
anisotropy  in  the  material  beneath  the  array. 


Chapter  4 

Velocities  Determinations  for  Propagation  to  Station  i 
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4.1  Introduction 

An  important  objective  of  this  research  was  to  use  the 
surface  waves  recorded  at  the  array  to  characterize  the 
geology  along  the  travel  path.  Most  of  the  path  from  the 
mines  to  the  array  lies  in  the  Animikie  basin  and  the  fold 
and  thrust  terrane  of  the  Penokean  orogen.  The  Keweenawan 
terrane  beneath  the  array  comprises  a  small  portion  of  the 
total  travel  path  (Figure  1.2).  The  propagation  to  the 
array  was  initially  modelled  by  Greenhalgh  (1979)  who  fit 
the  observed  travel  times  to  the  array  to  a  simple  two  layer 
model  (Figure  4.1). 

There  were  several  motives  for  studying  propagation  to 
the  array.  The  first  was  that  the  structure  of  the  Animikie 
basin  is  undetermined,  especially  compared  to  the  Keweenawan 
basin  which  lies  beneath  the  array.  The  basin  beneath  the 
array  has  been  studied  to  some  extent  by  several 
researchers.  Mooney  and  others  (1970)  performed  a 
refraction  investigation  which  defined  much  of  the  velocity 
structure  in  the  immediate  vicinity  of  the  array.  Mosher 
(1980)  studied  surface  wave  propagation  across  the  array  and 
this  research  has  expanded  upon  his  observations.  Chandler 
and  others  (1989)  interpreted  three  seismic  reflection 
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Figure  4 . 1 


Two  layer  model  obtained  by  Greenhalgh  to 
match  the  group  travel  times  from  the  mines 
to  Station  1 
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profiles  which,  although  not  in  the  exact  location  of  the 
array,  provide  the  most  detail  to  date  on  the  geologic 
structure  of  the  basins  which  flank  the  western  side  of  the 
Mid-Continent  rift.  The  structure  of  the  Animikie  basin,  on 
the  other  hand,  has  been  determined  in  less  detail. 

Recently,  the  geologic  map  of  the  Penokean  orogen  was 
revised  (Southwick  and  others,  1988) .  Previously,  the 
Animikie  basin  was  believed  to  pre-date  the  Penokean  orogen 
(Morey,  1983) .  The  current  interpretation  is  that  the 
Animikie  basin  is  a  foreland  basin  which  formed  as  a 
flexural  response  of  the  craton  to  loading  by  stacked  thrust 
sheets.  In  this  interpretation,  the  formation  of  the 
Animikie  basin  was  concurrent  with  the  tectonic  activity 
represented  by  the  fold  and  thrust  belt  (Southwick  and 
others,  1988) .  Although  the  basin  has  been  deformed,  the 
most  intense  deformation  occurred  in  the  internal  zone  of 
the  fold  and  thrust  belt  (Figure  1.5).  It  is  worthwhile, 
therefore,  to  try  to  determine  the  velocity  structure  of  the 
basin  and  interpret  it  in  terms  of  changes  in  thickness  of 
the  sedimentary  material  within  the  basin.  The  mine 
locations  provide  coverage  of  a  large  portion  of  the  basin. 

Another  motivation  for  focusing  the  study  on  this 
aspect  of  the  surface  wave  propagation  is  the  apparent 
lateral  refraction  and  scattering  which  occurs  for  these 
surface  waves.  The  easterly  azimuth  of  arrival  for  R2  makes 
it  evident  that  R2  has  been  either  reflected  from  the 
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Douglas  fault  or  scattered  by  a  short  wavelength 
heterogeneity  near  the  array.  Refraction  also  occurs  for  R1 
prior  to  arrival  at  the  array.  R1  arrives  first  at  Station 
1  for  all  source  areas.  For  most  of  the  mines,  this  can  be 
explained  by  the  refracting  effects  of  the  wedge  beneath  the 
array.  For  the  easternmost  mines,  Erie  and  Reserve,  a 
straight  path  to  the  array  includes  the  Keweenawan  wedge 
(Figure  4.2)  and  the  raypath  must  be  bent  for  R1  to  arrive 
at  Station  1  first.  Determining  the  amount  of  refraction 
which  occurs,  and  the  errors  in  computed  velocities  which 
result  from  the  refraction  is  important  for  characterizing 
the  suitability  of  these  types  of  data  for  other  seismic 
studies.  Again,  the  geometry  of  this  study,  with  several 
sources  spanning  an  azimuth  of  about  30  degrees  and  the 
seismic  waves  traversing  the  same  terrane,  is  appropriate 
for  this  objective. 

Station  1  was  singled  out  for  this  part  of  the  study. 
Since  it  is  located  over  basement,  the  Rl  arrivals  at 
Station  1  have  only  the  effects  of  propagation  through  the 
Penokean  terranes. 

4.2  Group  Velocities  to  Station  1 

The  best  approach  to  determine  the  effects  of 
propagation  through  the  Animikie  basin  is  to  first  assume 
the  simplest  model,  i.e.,  a  straight  path  from  source  to 
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Figure  4,2  Straight  raypaths  from  sources  to  Station  1. 

The  raypaths  for  Erie  and  Reserve  include  the 
Keweenawan  sedimentary  terrane. 
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receiver.  Inconsistencies  in  the  results  may  then  be 
explained  by  deviations  in  the  raypath  and  the  model  can  be 
revised  until  satisfactory  results  are  obtained.  Group 
velocities  to  Station  1  were  calculated  first.  The  group 
velocity  (U)  is  simply  the  distance  travelled  divided  by  the 
total  travel  time. 

U  =  x/t 

It  was  anticipated  that  the  group  velocities  determined  for 
Erie  and  Reserve  would  be  lower  than  the  velocities  for  the 
western  mines,  since  the  distance  used  would  be  too  low. 

To  determine  the  total  travel  time  of  the  surface 
waves,  the  blast  origin  times  had  to  be  estimated.  This 
could  be  accomplished  by  determining  the  P-wave  travel  time 
from  a  given  mine  to  the  array.  Greenhalgh  (1979)  measured 
the  travel  time  of  P-waves  to  each  station  of  the  array  for 
each  of  the  mines.  The  apparent  velocity  of  the  first  P- 
wave  arrival  at  Station  1  is  very  consistent  and  is  5.99 
km/s  for  all  mines  except  Erie,  which  has  an  apparent 
velocity  of  5.93  km/s.  Since  the  variations  in  apparent 
velocity  were  small,  a  reliable  estimate  of  the  origin  time 
could  be  made  by  measuring  the  first  break  on  Stations  l  and 
6.  Two  stations  were  used  for  quality  control.  The  first 
breaks  were  not  always  clear,  and  the  P-waves  are  very 
complicated  due  to  the  source  type  and  the  distance  from  the 
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source.  By  using  two  stations,  and  ensuring  that  the 
relative  time  of  the  first  motion  on  both  stations  fell 
within  a  prescribed  range,  the  estimate  of  the  blast  origin 
time  was  improved  by  a  few  seconds.  Station  1  was  used 
because  it  was  the  station  of  interest  and  Station  6  was 
used  because  the  seismograms  were  digitized  relative  to  the 
first  motion  on  that  station  so  it  was  part  of  the  velocity 
calculation.  The  actual  formula  used  to  calculate  group 
velocities  was: 

U  =  X/(Pi  +  +  t) 

where  P^  is  the  travel  time  of  the  P-wave  to 
Station  1 

Pj^_g  is  the  P-wave  travel  time  from  Station  l 
to  Station  6 

t  is  the  surface  wave  arrival  time  measured 
from  multiple  filter  analysis 

The  travel  times  were  averaged  for  each  mine  and  each  center 
frequency  between  0.4  Hz  and  1.5  Hz.  A  minimum  of  five 
readings  was  required  for  a  frequency  or  mine  to  be 
included.  Figure  4.3  shows  the  group  velocities  that  were 
determined.  The  results  show  excellent  agreement  between 
mines,  especially  in  light  of  the  ranges  observed  for  the 
array  group  velocities  (Figure  4.4).  The  number  of  records 
used  for  each  mine  was  between  six  and  thirteen.  The  most 
consistent  travel  times  for  each  mine  were  at  0.8  Hz  and  1 


GROUP  VELOCITIES  MEASURED  AT  STATION  1 
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Figure  4.3  Group  velocity  of  R1  as  a  function  of  frequency,  mea.sured  at  Station 
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Figure  4.4  Comparison  of  R1  group  velocities  to  Station  1  and  across  the  CMSA 
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Hz.  The  standard  deviations  for  the  computed  velocities  at 
these  frequencies  were  less  than  0.1  km/s.  The  range  of 
velocities  is  from  about  2.65  km/s  to  2.85  km/s.  Erie  and 
Reserve  are  slower  than  the  mines  to  the  west,  but  probably 
not  significantly  so.  Butler  should  have  the  highest 
velocity  since  its  actual  raypath  should  have  the  least 
amount  of  deviation  from  its  assumed  raypath.  This  is  not 
the  case.  U.S.  Steel  has  the  highest  group  velocity.  This 
may  indicate  a  localized  velocity  high  associated  with  that 
source  area. 

The  next  step  in  the  group  velocity  analysis  was  to 
separate  the  velocities  into  an  Animikie  basin  velocity  and 
a  fold  and  thrust  terrane  velocity.  Since  the  fold  and 
thrust  terrane  is  extremely  complicated,  but  of  a  high 
metamorphic  grade,  it  was  assigned  a  constant  velocity  of 
3.25  km/s.  The  velocity  for  the  Animikie  basin  may  be 
calculated  by  determining  the  percentage  of  the  raypath  in 
each  terrane.  The  group  velocity  for  the  total  path  is  a 
simple  composite  of  the  travel  time  through  each  terrane: 

Utotal  =  (^1  +  X2)/(Xi/Ui  +  X2/U2)  (Knopoff,1969) 

where  1  and  2  refer  to  the  respective  terranes. 
Interestingly,  the  path  length  through  the  fold  and  thrust 
terrane  is  approximately  the  same  length  for  all  of  the 
mines.  It  varies  between  76  and  80  kilometers.  Because  the 
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path  length  only  varies  through  the  Animikie  basin,  no 
changes  in  the  relative  position  of  the  different  source 
areas  on  the  velocity  curve  occurs.  Figure  4.5  shows  the 
velocities  derived  for  the  basin  using  this  approximation. 
The  velocities  range  between  2.2  km/s  and  2.6  km/s.  The 
only  significant  change  in  the  curve  as  compared  to  the 
total  velocity  curve  is  that  the  high  velocity  associated 
with  U.S.  Steel  is  more  pronounced.  The  change  corresponds 
to  about  6  seconds  difference  in  total  travel  time.  This  is 
probably  close  to  the  limit  of  detectability  for  a  velocity 
anomaly.  The  standard  deviation  for  U.S.  Steel  group 
velocities  was  around  .06  km/s.  If  the  errors  in  velocity 
are  twice  the  standard  deviation,  anything  larger  than  .12 
km/s  is  significant.  From  this  perspective,  at  the  source 
areas  adjacent  to  U.S.  Steel,  the  increase  in  velocity  at  1 
Hz  is  insignificant  between  Hibbing  and  U.S.  Steel,  since  it 
is  exactly  .12  km/s.  It  is  significant  between  U.S.  Steel 
and  Erie,  since  the  difference  is  .21  km/s. 

4.3  Two  Source  Phase  Velocities 

The  correlation  coefficients  are  very  high  between  R1 
arrivals  recorded  on  Station  1  from  different  mines.  This 
can  be  exploited  to  obtain  velocity  information  by  making 
assumptions  about  what  gives  rise  to  this  high  correlation. 

A  high  correlation  coefficient  implies  that  the  phase  shift 


GROUP  VELOCITIES  ESTIMATED  FOR 
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Figure  4.5  R1  group  velocity  as  a  function  of  frequency  estimated  for  the  Animikie 
Basin 
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between  the  two  records  being  compared  is  a  pure  delay. 

This  implies  that  the  waves  have  either  propagated  through 
identical  media,  or  taken  very  similar  paths  from  source  to 
receiver.  It  was  decided  that  the  records  which  had  very 
high  correlation  coefficients  could  be  used  to  examine  these 
two  possibilities  for  the  propagation  of  Rl.  From  the 
cross-correlation  coefficients,  phase  velocities  could  be 
computed  using  two  sources  and  taking  the  difference  in 
travel  time  of  a  given  phase  to  determine  the  near  source 
phase  velocity. 

The  two  source  method  is  somewhat  unorthodox,  since 
phase  velocities  are  usually  measured  at  two  stations  along 
the  great  circle  path  from  a  single  source  (Bloch  and  Hales, 
1968) .  By  using  two  sources,  there  are  implicit  assumptions 
that  are  necessary  to  make  the  measurement  valid.  The  first 
is  that  the  material  being  sampled  along  the  two  raypaths 
has  the  same  dispersion.  The  second  assumption  is  that  the 
measured  phase  velocity,  which  is  assigned  to  the  source 
area,  is  constant  for  all  frequencies. 

The  phase  shifts  were  measured  by  multiple  filter 
analysis.  Prior  to  filtering,  the  maximum  value  of  the 
cross-correlation  function  was  determined  for  the  two 
records  being  compared.  The  cross-correlation  function  had 
to  attain  a  maximum  value  of  at  least  0.8  for  two  records  to 
be  used.  Many  records  had  cross-correlation  coefficients 
above  0.9.  In  addition,  the  same  conditions  on  the  estimate 
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of  origin  time  that  were  used  to  obtain  group  velocities  to 
the  array  were  imposed.  These  two  conditions  limited  the 
number  of  records  available.  Most  of  the  mines  had  five  or 
more  records  which  met  these  conditions.  Phase  shifts  were 
measured  from  0.4  Hz  to  2.0  Hz  at  0.2  Hz  increments  between 
center  frequencies.  Figure  4.6  shows  an  envelope  derived 
from  the  multiple  filter  analysis.  It  is  apparent  that  the 
method  does  not  work  for  frequencies  above  1.2  Hz. 

The  first  result  obtained  was  that  the  condition  that 
the  correlation  coefficient  attain  a  value  of  0.8  eliminated 
all  mines  east  of  U.S.  Steel.  No  correlations  above  0.8 
were  obtained  for  records  from  Erie  and  Reserve.  This  is 
probably  due  to  the  heterogeneous  nature  of  the  eastern 
Animikie  basin. 

The  phase  differences  between  the  western  mines  were 
obtained  and  averaged  for  the  frequency  range  0.4  to  1.2  Hz. 
To  convert  to  velocity,  the  raypath  must  be  assumed.  Two 
end  points  were  tested.  The  first  assumes  a  straight 
raypath  from  source  to  receiver.  The  phase  velocities  are 
shown  in  Figure  4.7.  Some  interesting  features  emerge.  The 
curves  are  generally  flat,  and  therefore  fit  the  assumption 
of  constant  phase  velocity  relative  to  frequency  for  most 
mine  pairs.  It  is  interesting  that  the  Butler-National 
curve  and  the  National-Hibbing  curve  are  not  flat  but  the 
Butler-Hibbing  curve  averages  the  two  and  is  very  flat. 

There  is  very  good  evidence  for  a  high  velocity  associated 
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Figure  4.7  Two  source  phase  velocity  as  a  function  of  frequency  assuming  a  straight 
raypath  for  R1  to  Station  1 
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with  U.S.  Steel.  The  local  velocity  between  U.S.  Steel  and 
Hibbing  is  above  3  km/s.  The  anomaly  drops  to  about  2.75 
km/s  when  comparing  Butler  and  U.S.  Steel.  Butler  and  U.S. 
Steel  are  about  26  kilometers  apart,  so  this  is  consistent 
with  the  high  velocity  being  localized  near  U.S.  Steel. 

A  possible  source  for  the  high  velocities  observed  near 
U.S.  Steel  is  a  high  velocity  material  within  the  basement 
near  the  mine.  The  Animikie  Group  is  thin  within  the  Mesabi 
Range,  so  the  near  source  velocities  are  mostly  sensitive  to 
basement  lithology.  Within  the  northeastern  Animikie  basin, 
a  broad  gravity  high  exists  and  has  been  interpreted  as  a 
greenstone  belt  within  the  Archean  basement.  The  maximum 
values  for  the  anomaly  occur  just  south  of  U.S.  Steel,  very 
near  the  source  area  (Chandler,  1985) .  The  greenstone  belt 
contains  metavolcanics  which  are  likely  to  have  higher 
velocities  than  the  felsic  rocks  of  the  Giants  Range 
batholith  which  underlies  much  of  the  northern  Animikie 
basin. 

A  possibility  that  was  considered  at  the  time  these 
phase  velocities  were  calculated  was  that  the  Rl  surface 
waves  were  skirting  the  edge  of  the  basin  and  arriving  at 
the  array  from  a  northwesterly  azimuth  (Figure  4.8).  This 
was  based  on  the  lack  of  evidence  for  refraction  increasing 
across  the  array  for  Rl  from  the  array  phase  velocities. 

This  type  of  raypath  would  explain  completely  the  high 
correlation  coefficients  for  Rl  at  Station  1  since  Rl  would 


Figure  4 . 8  Raypath  used  to  calculate  near  source  phase 
velocities.  R1  travels  along  strike  of  the 
Mesabi  Range  as  it  leaves  the  source  area. 
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have  a  nearly  identical  raypath  for  all  source  areas.  To 
test  this,  t.ie  phase  velocities  were  calculated  assuming  the 
rays  travelled  along  the  strike  of  the  Mesabi  Range.  These 
phase  velocities  are  shown  in  Figure  4.9. 

The  phase  velocities  quickly  take  on  unrealistically 
high  values  for  this  case.  Therefore,  the  straight  path 
approximation  is  closer  to  correct.  The  narrow  range  of 
group  velocities  measured  for  this  case  also  supports  this. 
However,  the  arrival  of  R1  at  Station  1  has  still  not  been 
satisfactorily  explained  by  the  straight  path  assumption. 

The  raypaths  may  converge  within  the  Animikie  basin.  This 
is  due  to  lateral  velocity  gradients  present  within  the 
basin.  The  deviation  from  the  straight  raypath  is  not 
extreme,  but  should  increase  for  source  areas  to  the  east. 
The  near  source  high  velocity  associated  with  U.S.  Steel  may 
mask  this  effect  for  that  mine.  The  raypaths  from  Erie  and 
Reserve  should  show  the  most  deviation  from  a  straight  path. 
The  results  from  modelling  raypaths  and  travel  times  are 
discussed  in  Chapter  5. 

4 . 4  Summary  of  Results  from  Velocity  Analysis  at  Station  1 

The  velocities  from  the  mines  to  the  array  were 
measured  for  various  raypaths  to  Station  1.  The  group 
velocities  to  Station  1  assumed  a  straight  raypath  from  each 
mine.  The  group  velocities  obtained  are  very  consistent  and 
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range  from  2.64  km/s  to  2.85  ]cm/s.  This  range  is  lower  than 
the  range  of  group  velocities  determined  for  the  array 
(Figure  4.4).  There  is  evidence  of  higher  than  average 
group  velocities  for  records  from  U.S.  Steel.  Group 
velocities  were  estimated  for  the  Animikie  basin  and  range 
from  2.2  km/s  to  2.6  km/s. 

Near  source  phase  velocities  were  determined  for 
Station  1  records  which  had  correlation  coefficients  above 
0.8.  The  phase  velocities  determined  by  assuming  that  R1 
follows  a  straight  raypath  to  Station  1  range  from  1.8  km/s 
to  3.6  km/s.  The  phase  velocity  measurements  also  confirm 
the  presence  of  a  velocity  high  near  U.S.  Steel.  The  phase 
velocities  were  also  calculated  assuming  a  raypath  which 
skirts  the  northern  edge  of  the  Animikie  basin.  These 
velocities  range  from  3.10  km/s  to  7.07  km/s.  The  actual 
raypath  is  somewhere  between  these  two  ranges,  but  is  closer 
to  the  straight  raypath. 


Chapter  5 

Two  Dimensional  Ravtracina 
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5.1  Introduction 

In  an  effort  to  further  understand  the  propagation  of 
the  R1  arrival  to  the  array,  two  dimensional  models  were 
constructed  to  fit  velocities  and  raypaths  for  R1  through 
the  major  Precambrian  tectonic  terranes.  The  velocity 
structure  was  designed  to  match  the  group  travel  times  to 
the  array  stations  for  a  1  Hz  arrival.  The  travel  time 
observations  are  most  consistent  at  0.6  Hz,  0.8  Hz  and  1.0 
Hz.  The  travel  times  for  these  frequencies  differ  only  by  a 
few  seconds  between  each  frequency,  so  an  initial  model  for 
a  single  frequency  is  adequate  to  constrain  which  features 
most  affect  the  R1  arrival. 

There  were  several  constraints  on  the  models  arising 
from  the  velocity  analyses  presented  in  previous  chapters. 
There  were  also  several  measurements  of  near  surface  P-wave 
velocities  available  from  refraction  studies  (Mooney  and 
others,  1970) .  A  rough  estimate  of  the  Rayleigh  wave 
velocity  can  be  obtained  from  the  P-wave  velocity  by 
calculating  the  velocity  of  a  Rayleigh  wave  in  a  uniform 
elastic  medium: 

C  *  .9194  Vp/(3)^/^  (Ewing  and  others,  1957) 
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The  available  constraints  on  the  model  are  as  follows; 


(1)  The  direct  path  group  velocity  to  Station  1  is 
about  2.75  km/s. 

(2)  The  average  group  velocity  for  the  Animikie  basin 
has  been  estimated  as  2. 2 -2. 6  km/s. 

(3)  The  phase  velocity  near  the  source  has  been 
estimated  as  approximately  2.8  km/s  for  the 
western  source  areas.  A  local  velocity  high  is 
associated  with  U.S.  Steel.  The  velocity  near 
U.S.  Steel  may  be  above  3.0  km/s.  No  near  source 
velocities  were  obtained  for  Erie  and  Reserve. 

(4)  The  R1  group  velocity  across  the  array  is  1.7 
km/s.  The  measured  azimuth  for  R1  is  320  degrees. 
It  is  also  known  however  that  the  R1  arrival  is 
always  observed  on  Station  1  first,  and  therefore 
does  not  follow  a  straight-line  path  from  the 
mines  to  the  other  array  receivers. 

(5)  The  P-wave  velocities  which  have  been  measured 
are; 


5.0  -  5.5  km/s  for  Keweenawan  volcanics 
5.7  -  6.2  km/s  for  pre-Keweenawan  felsic  to 
intermediate  intrusive  rocks,  ie,  Pcnokean 
intrusives 

6.6  -  7.1  km/s  for  Keweenawan  mafic  intrusives, 

ie,  Duluth  Complex  rocks 

2.74  km/s  for  the  Hinckley  Sandstone 

3.65  km/s  for  the  upper  Fond  du  Lac  Formation 
(Mooney  and  others,  1970,  Greenhalgh,  1979, 
Chandler  and  others,  1989) 

These  values  yield  the  following  initial  estimates 
for  Rayleigh  wave  velocities; 

2.65  -  3.0  km/s  for  Keweenawan  volcanics 

3.0  -  3.3  km/s  for  Penokean  crystalline  rocks 

3.5  -  3.7  km/s  for  Duluth  Complex  rocks 

1.4  -  1.9  km/s  for  western  basin  sedimentary  rocks 
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5.2  The  RAYS 4  Program 

Since  surface  waves  travel  along  the  free  surface,  it 
is  possible  to  consider  the  terranes  as  two  dimensionally 
varying  velocity  structures.  Cerveny  and  others  (1977) 
describes  a  forward  modelling  approach  which  can  be  applied 
to  surface  waves  in  laterally  varying  media.  For  a  series 
of  surface  locations,  the  local  phase  velocity  may  be 
determined  if  the  variation  of  the  elastic  parameters  X,  n, 
p  are  known  as  a  function  of  depth.  If  these  are  known, 
local  phase  velocities  can  be  determined  at  discrete 
frequencies.  For  each  individual  frequency,  standard  two 
dimensional  raytracing  methods  may  be  used  to  find  the 
raypath  for  the  individual  frequencies.  The  only 
measurements  available  for  the  modelling  are  the  group 
travel  times,  from  which  phase  travel  times  can  be 
approximated,  and  some  regional  measurements  of  density  and 
P-wave  velocities.  The  raypath  is  constrained  in  that  Rl 
must  arrive  at  the  array  and  it  must  arrive  at  Station  1 
first.  If  the  travel  times  are  matched,  the  resulting  two 
dimensional  velocity  structure  can  then  be  interpreted  in 
terms  of  changes  in  thickness  of  material  or  changes  in 
lithology.  A  complete  analysis  would  also  include  synthesis 
of  three  dimensional  geologic  structure  and  the  resulting 
dispersion  curves  for  Rl.  However,  this  is  very  time  and 
computer  intensive,  and  the  quality  of  this  data  set  does 
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not  merit  this  rigorous  treatment.  Also,  the  geologic 
constraints  are  not  available  to  make  this  a  realistically 
tractable  problem. 

The  raytracing  was  done  using  a  program  called  RAY84, 
an  interactive  two  dimensional  raytracing  program  designed 
for  use  with  body  waves.  The  travel  time  and  position  in 
time  for  the  rays  are  computed  from  a  system  of  first  order 
differential  equations: 

dX(t)/dt  »  V(X,Z)  sin(e) 

dZ(t)/dt  =  V(X,Z)  cos(e) 

d0(t)/dt  =  (dV/dX)  cos (8)  -  (dV/dZ)  sin (6) 

where  0  is  the  initial  angle  of  the  ray.  The  velocity  model 
V(X,Z)  is  supplied  as  a  data  file  and  the  initial  values  for 
X,  Z,  0,  and  t  are  supplied  interactively.  The  program 
simultaneously  integrates  the  three  equations  over  small 
steps  in  time  (Luetgart,  1988) . 

The  modelling  was  done  on  a  Microvax/VMS  workstation  at 
the  University  of  Wisconsin  Oshkosh.  Access  to  the 
facilities  was  provided  by  Dr.  John  Karl.  The  program  was 
adapted  for  this  purpose  by  treating  the  source  locations, 
or  north,  as  the  negative  Z  direction  and  south  as  the 
positive  Z  direction.  The  positive  X  direction  was  east. 

The  model  boundaries  are  defined  by  interfaces  extending 
from  west  to  east.  These  represent  lithologic  interfaces. 
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These  layers  are  allowed  to  pinch  out.  The  only  difficulty 
arose  in  defining  vertical,  or  north-south,  contacts,  such 
as  the  eastern  boundary  of  the  Animikie  basin.  This  was 
overcome  by  taking  advantage  of  the  manner  in  which  the 
velocities  may  be  defined  within  the  model.  Velocities  are 
defined  along  vertical  grid  lines  which  cross  all  lithologic 
boundaries.  Wherever  a  velocity  grid  line  crosses  a 
lithologic  interface,  the  velocity  above  and  below  the 
interface  is  defined.  This  allows  lateral  and  vertical 
velocity  gradients  to  be  used  in  the  model  as  well  as  first 
order  vertical  velocity  discontinuities  (Luetgart,  1988) . 

By  placing  two  closely  spaced  velocity  grid  lines  near  the 
contact  between  the  Duluth  Complex  and  the  Animikie  basin 
and  assigning  a  velocity  of  3.5  km/s  for  the  easternmost 
grid  line,  the  contact  could  be  simulated  by  a  steep  lateral 
velocity  gradient  within  the  basin.  The  Duluth  Complex  was 
incorporated  into  the  model  not  as  a  lithologic  boundary  but 
by  a  high  lateral  velocity  gradient  between  itself  and  the 
Animikie  basin. 

Figure  5.1  shows  the  lithologic  boundaries  and  the 
locations  of  the  velocity  grid  lines  used  to  model  Rl 
propagation.  Several  models  are  presented  in  the  next 
section  along  with  the  important  features  that  most  strongly 
influence  the  raypaths  to  the  array.  Many  of  the  short 
wavelength  spatial  features  from  the  geologic  map  have  been 
smoothed  since  they  lead  to  reflection  and  scattering  of  the 
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Lithologic  boundaries  and  velocity  grid  lines 
used  to  model  the  Precambrian  geologic 
terranes  with  the  RAYS 4  program.  Source  and 
receiver  locations  are  also  indicated 


Figure  5 . 1 
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waves.  The  multipathing  of  R1  is  more  likely  caused  by 
large  scale,  smoothly  varying  velocity  changes  within  a 
terrane,  or  by  refraction  upon  crossing  the  boundary  between 
two  tectonic  terranes.  The  three  tectonic  terranes  which 
affect  the  Rl  propagation  are;  the  Animikie  basin,  the  fold 
and  thrust  terrane,  and  the  Keweenawan  basin  which  flanks 
the  St.  Croix  Horst  to  the  west. 


5.3  Results  from  Ravtracina 

Three  models  are  presented.  The  first  two  models 
evolved  into  the  final  model  so  presenting  them  helps 
clarify  features  of  the  final  model. 

The  first  model  attempted  to  determine  if  Rl  arrives  at 
the  array  simply  by  refracting  where  it  crosses  terrane 
boundaries.  For  this  case,  each  terrane  was  assigned  a 
constant  velocity.  The  Animikie  basin  was  assigned  a 
velocity  of  2.5  km/s,  the  fold  and  thrust  terrane,  3.15 
km/s,  and  the  Keweenawan  basin,  1.7  km/s. 

A  feature  of  the  RAY84  program  which  was  used 
extensively  in  the  modelling  process  Is  che  option  to  set  a 
target  in  the  "subsurface”.  The  target  is  defined  by  X-Z 
coordinates  and  radius.  All  rays  which  encounter  the  target 
are  terminated  upon  encountering  it.  By  locating  the  target 
at  Station  6  of  the  array,  and  giving  the  target  a  radius  of 
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14  kilometers,  all  rays  which  cross  the  aperture  of  the 
array  can  be  found. 

Three  important  results  were  obtained  from  the  constant 
velocity  model,  two  of  which  were  anticipated.  The  first  is 
that  the  refracting  effect  of  the  basin  beneath  the  array  is 
considerable  (Figure  5.2).  The  angle  at  which  Rl  crosses 
the  array  is  practically  independent  of  the  angle  of 
incidence  at  the  boundary  because  of  the  high  velocity 
contrast  between  the  Penokean  age  rocKS  and  the  Keweenawan 
rocks . 

The  second  result  is  that  with  no  lateral  velocity 
changes  within  the  Animikie  basin,  Rl  does  not  reach  the 
array  from  eastern  source  areas  such  as  Erie  (Figure  5.3) . 
This  is  due  to  both  the  refracting  effect  of  the  western 
basin,  and  the  geometry  of  the  Mid-Continent  rift  system. 

The  low  velocities  in  the  basin  prevent  any  direct  path 
surface  waves  from  reaching  the  array  and  since  the  strike 
of  the  rift  system  changes  to  a  more  southerly  orientation 
in  the  vicinity  of  the  array,  all  rays  which  leave  the 
eastern  source  areas  miss  the  array.  Therefore  a  lateral 
gradient  is  required  in  the  Animikie  basin  which  will  bend 
the  raypaths  to  the  south  as  they  near  the  southern  boundary 
of  the  basin. 

The  third  effect,  which  was  not  anticipated,  is  the 
strong  effect  that  the  shape  of  the  southern  boundary  of  the 
Animikie  basin  has  on  the  raypaths.  The  shape  of  this 
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Figure  5.2  Raytracing  results  for  model  with  constant  velocity  in  each  terrane.  The  source  is 
Hibbing  The  refracting  effect  of  the  Keweenawan  basin  can  be  seen 


Figure  5.3  Raytracing  results  for  model  with  constant  velocity  in  each  terrane 
The  source  is  Erie.  No  rays  encounter  the  array 
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boundary,  as  modelled  from  the  geologic  map,  creates  a 
shadow  zone  around  the  array  for  rays  originating  from 
Butler  and  National  (Figure  5.4).  This  effect  could  not  be 
compensated  for  without  changing  the  direction  of  concavity 
of  the  boundary.  If  R1  propagates  as  a  wavefront  from  the 
mines  to  the  array,  it  must  cross  this  boundary  without  the 
divergence  in  the  raypath  that  this  model  indicates. 
Unfortunately,  this  part  of  the  model  is  the  most  poorly 
constrained.  It  is  known  that  the  southern  part  of  the 
basin  has  undergone  more  deformation  than  the  northern  part 
(Morey,  1983,  Holst,  1984,  Southwick  and  others,  1988),  but 
it  is  also  reasonable  to  assume  that  since  the  metamorphic 
grade  is  higher  in  the  fold  and  thrust  terrane,  the  velocity 
increases  in  that  terrane.  At  this  time,  it  is  not  known 
what  the  order  of  the  velocity  changes  are  between  the 
Animikie  basin  and  the  fold  and  thrust  terrane,  or  where 
they  occur. 

It  is  possible  that  the  transition  from  the  Animikie 
basin  to  the  fold  and  thrust  terrane  affects  the  surface 
waves  significantly.  The  fold  and  thrust  terrane  is 
geologically  complicated  and  could  contain  short  wavelength 
features  which  may  act  as  scatterers.  Rl  may  be  regenerated 
within  this  terrane.  The  distance  from  the  northern 
boundary  of  the  fold  and  thrust  terrane  to  the  array  is 
about  80  kilometers.  Mosher  (1980)  determined  the  distance 
at  which  wavefront  curvature  could  be  detected  by  the  array 
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Figure  5.4  Raytracing  results  for  model  with  constant  velocity  in  each  terrane. 

The  source  is  National.  The  shape  of  the  southern  boundary  of  the  Animikie 
Basin  creates  a  shadow  zone  which  causes  most  rays  to  miss  the  array 
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to  be  3  array  radii,  or  approximately  42  kilometers. 
Therefore  it  is  not  possible  to  determine  if  Rl  is  scattered 
by  some  feature  within  the  fold  and  thrust  terrane.  This 
prospect  could  be  further  studied  through  normal  mode 
synthesis  of  an  Rg  wave  and  two  or  three  dimensional 
modelling  of  its  transition  between  terranes,  but  more 
knowledge  of  the  P  and  S  wave  velocities  and  densities  of 
the  materials  is  required. 

A  second  model  was  based  on  the  structure  of  the 
Animikie  basin  obtained  by  Ferderer  (1988)  (Figure  5.5). 

The  contours  in  Figure  5.5  represent  possible  thicknesses  of 
the  basin  obtained  from  Werner  deconvolution  of 
aeromagnetic  data.  The  general  structure  has  a  gently 
sloping  shelf  in  the  northern  part  of  the  basin,  and 
considerable  thicknesses  of  material  present  to  the 
southeast.  According  to  Ferderer  (1988)  the  dips  obtained 
for  this  model  may  correspond  to  dips  within  the  underlying 
basement,  or  they  may  represent  the  dip  of  the  iron 
formation  within  the  Animikie  basin  which  gives  rise  to  the 
magnetic  signal.  A  third  possibility  discussed  by  Ferderer 
for  his  dips  is  a  strong  remanent  magnetization  (rather  than 
induced  magnetization)  in  the  unit  which  is  the  signal 
source . 

The  second  model  contoured  the  velocity  structure  of 
the  basin  based  on  the  depth  contours  in  Figure  5.5.  Each 
contour  represents  a  second  order  decrease  in  velocity.  The 
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northern  boundary  of  the  basin  was  assigned  a  velocity  of  3 
km/s.  The  velocity  decreased  to  a  value  of  2  km/s  at  the 
6000  meter  contour  of  Figure  5.5.  The  velocities  were 
estimated  from  the  ranges  of  measured  velocities.  The  fold 
and  thrust  terrane  was  given  a  velocity  of  3.25  km/s.  The 
Keweenawan  basin  was  given  a  velocity  which  changed  from  1.8 
km/s  at  the  western  edge  to  1.4  km/s  at  the  Douglas  Fault. 

The  first  source  area  tried  was  Hibbing  (Figure  5.6). 
This  model  matched  the  arrivals  from  Hibbing  to  within  1 
second  of  the  measured  travel  times.  When  National  and 
Butler  were  used  as  source  areas,  the  divergence  of  the  rays 
was  even  greater  than  for  the  constant  velocity  model 
(Figure  5.7).  Another  problem  with  this  mocel  was  that  the 
southern  array  stations  were  difficult  to  hit.  The  southern 
boundary  of  the  Animikie  basin  was  smoothed  and  the  4000  and 
5000  meter  contours  were  changed  to  mimic  the  3000  meter 
contour  in  an  attempt  to  eliminate  the  shadow  zone  (Figure 
5.8).  When  this  was  done,  a  new  shadow  zone  was  created  for 
rays  originating  at  Hibbing.  Eventually,  it  was  decided 
that  this  model  was  too  complicated  and  it  was  abandoned  in 
favor  of  a  simpler  model.  The  results  from  Hibbing  did 
yield  good  estimates  for  the  velocity  contrasts  needed  to 
match  the  travel  times  to  the  array,  and  the  final  model 
incorporates  the  overall  trends  obtained  by  Ferderer  (1988) 
for  the  Animikie  basin. 
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Figure  5.7  Raypaths  for  model  based  on  Ferderer  (1988).  For  rays  originating  at 

Butler  and  National,  the  shadow  zone  around  the  array  is  even  larger  than 
for  the  constant  velocity  model. 
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Figure  5.8  Raypaths  for  National  obtained  by  smoothing  the  contours  of  the  previous 
model.  The  northern  array  stations  are  reached. 
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The  third  model  used  a  velocity  structure  for  the 
Animikie  basin  in  which  velocity  generally  decreases  to  the 
south  and  east.  The  final  model  is  shown  in  Figure  5.9  and 
the  resulting  raypaths  are  shown  in  Figure  5.10.  The 
southern  boundary  of  the  Animikie  basin  was  smoothed 
considerably.  From  a  geophysical  aspect,  this  may  be  more 
plausible  than  modelling  the  actual  geologic  boundary. 

Since  the  terranes  are  tectonically  related,  the  velocity 
variations  may  be  more  smooth  than  the  geologic  contacts. 

In  the  final  model,  the  velocities  within  the  basin  were 
defined  only  at  the  northern  and  southern  boundaries  and 
lateral  changes  were  incorporated  through  placement  of  the 
velocity  grid  lines.  The  minimum  velocity  used  for  the 
basin  was  1.9  km/s.  The  fold  and  thrust  terrane  was 
initially  assigned  a  constant  velocity  of  3.0  km/s,  but  the 
rays  were  arriving  2  to  4  seconds  early  at  all  of  the  array 
stations.  A  negative  velocity  gradient  from  northwest  to 
southeast  with  a  velocity  range  from  3.15  km/s  at  the 
northern  boundary  to  2.9  km/s  at  the  contact  between  the 
Penokean  crystalline  material  and  the  western  Keweenawan 
basin  was  introduced  into  the  fold  and  thrust  terrane  in 
order  to  delay  the  travel  times  and  bend  the  raypaths  to  a 
more  easterly  azimuth.  The  velocity  was  lowered  in  the  fold 
and  thrust  terrane  terrane  because  the  azimuth  across  the 
array  was  330  degrees  for  the  model,  and  the  observed 
azimuth  is  320  degrees.  This  velocity  structure  for  the 


Final  two-dimensional  model  for  1  Hz  R1  surface  wave 
obtained  from  raytracing 
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fold  and  thrust  terrane  can  be  justified  geologically  by 
associating  the  decrease  in  velocity  with  the  presence  of 
more  felsic  material  in  the  eastern  parts  of  the  fold  and 
thrust  terrane. 

Although  the  lower  velocity  within  the  fold  and  thrust 
terrane  corrected  the  travel  times,  the  azimuth  appears  to 
be  solely  controlled  by  the  velocity  of  the  Keweenawan 
basin.  Unless  a  gradient  is  incorporated  into  the  basin, 
the  rays  cross  the  basin  at  a  constant  azimuth  of  330 
degrees,  independent  of  the  angle  of  incidence.  In  the 
example  in  Fxgure  5.8,  in  which  the  velocity  varies  from  1.8 
km/s  to  1.4  km/s,  the  rays  enter  the  basin  with  an  azimuth 
of  330  degrees  and  encounter  the  Douglas  Fault  at  an  azimuth 
of  309  degrees. 

Since  the  measured  azimuth  using  all  six  array  stations 
and  various  station  combinations  is  lower  than  can  be 
obtained  by  using  a  constant  velocity  for  the  Keweenawan 
basin,  there  must  be  a  velocity  gradient  within  the 
Keweenawan  basin  with  velocity  decreasing  as  the  wedge 
thickens.  Although  a  change  in  azimuth  across  the  array  was 
only  observed  for  group  velocities,  it  was  predicted  and  is 
necessarily  present  to  obtain  the  correct  six  station 
azimuth  for  the  raytracing  model. 
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5 . 4  Two  Laver  Dispersion  Models  for  the  Animikie  Basin 

After  the  raytracing  was  completed,  new  dispersion 
curves  were  generated  for  the  average  group  velocity  of  the 
Animikie  basin  as  described  in  Chapter  i .  An  average 
velocity  of  3.0  km/s  was  used  for  the  fold  and  thrust 
terrane  and  straight  raypaths  through  the  basin  were  assumed 
since  these  are  approximately  correct  as  shown  by  the  final 
model  and  the  results  of  Chapter  4.  The  resulting  curves 
are  shown  in  Figure  5.11. 

Next,  theoretical  group  velocity  versus  period  curves 
were  generated  for  a  two  dimensional  layered  structure  in  an 
attempt  to  determine  the  average  thickness  and  velocities  of 
the  Animikie  basin.  The  curves  were  computed  from  a 
synthetic  seismogram  program  written  by  Wang  (1985  a,  b,  c, 
d) .  The  program  computes  phase  and  group  velocity  curves 
for  any  number  of  layers  over  an  elastic  half  space  for  as 
many  modes  as  requested.  This  particular  effort  only 
attempted  to  model  two  layers,  the  Animikie  Group  and  an 
underlying  half  space  which  represents  basement.  The  layer 
representing  the  Animikie  Group  may  be  subsequently 
thickened  or  thinned  and  the  group  velocity  at  1  Hz 
determined  in  an  attempt  to  match  the  trends  modelled  by  the 
raytracing. 

Greenhalgh  (1979)  made  a  two  dimensional  model  for  the 
entire  structure  between  the  mines  and  Station  1  (Figure 
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Figure  5.11  Estimate  of  the  Animikie  Basin  group  velocity  vs  period  based  on  the 

assumptions  of  straight  raypaths  through  the  basin  and  a  lower  velocity  in 
the  fold  and  thrust  terrane  than  was  used  in  Chapter  4 . 
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4.1).  P-wave  velocities  were  obtained  from  refraction 
studies  and  are  also  reported  in  Greenhalgh  (1979).  For  the 
Animikie  Group  the  P-wave  velocity  is  between  4.5  and  5.0 
km/s.  The  underlying  rocks  have  a  P-wave  velocity  between 
5.5  and  6.0  km/s.  These  were  the  only  constraints  available 
for  the  model  by  Greenhalgh.  The  S-wave  velocities  used  by 
Greenhalgh  were  obtained  by  assuming  Poisson's  ratio  equal 
to  0.25  for  both  layers.  The  densities  he  used  were  2.32 
g/cm^  for  the  upper  layer  and  2.64  g/cm^  for  the  half  space. 

The  curve  for  the  Animikie  basin  could  not  be  matched 
very  well  using  Greenhalgh 's  model.  Many  models  were 
attempted  within  the  range  of  possible  velocities.  The  best 
fit  was  obtained  by  raising  the  densities  in  both  layers  and 
using  Greenhalgh 's  velocities  and  a  layer  thickness  of  1.2 
kilometers  (Figure  5.12).  The  densities  used  were  taken 
from  Carlson  (1985).  The  best  fit  was  obtained  with  a 
density  for  the  Animikie  Group  of  3.1  g/cm^.  Carlson 
reported  a  density  of  3.45  g/cm^  for  the  Biwabik  Formation 
which  is  an  iron  formation  within  the  Animikie  Group.  She 
also  reported  a  density  of  3.07  g/cm^  for  the  non-magnetic 
rocks  of  the  Cuyuna  Range,  which  though  older,  should  have 
densities  similar  to  the  Animikie  Group  since  the 
lithologies  are  similar. 

The  density  used  for  the  Animikie  Group  may  be  somewhat 
high  (Chandler,  personal  communication).  If  the  data  at  0.4 
Hz  are  not  considered,  a  better  fit  could  be  obtained  with  a 
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Figure  5,12  Comparison  of  theoretical  group  velocity  curves  for  a  single  layer  over 
half  space  with  measured  values.  These  curves  illustrate  the  effect  of 
Increasing  the  density  in  the  layer. 
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lower  density  in  the  upper  layer  (Figure  5.12).  The  0.4  Hz 
velocity  values  had  large  standard  deviations,  and  may  not 
be  as  reliable  as  the  higher  frequency  values. 

The  basin  is  underlain  to  the  north  by  the  Giants  Range 
Batholith  for  which  Carlson  reported  densities  between  2.59 
and  2.84  g/cm^.  There  may  also  be  some  Archean  greenstones 
underlying  the  basin  with  densities  between  2.65  g/cm^  and 
3.14  g/cm^.  The  best  dispersion  curve  was  obtained  with  a 
density  in  the  half space  of  2.7  g/cm^.  This  is  in 
agreement  with  either  rock  type  being  present. 

This  two  layer  model  fits  the  data  very  well.  It  also 
agrees  with  Chandler's  (1983)  model  of  the  northern  Animikie 
basin  in  which  a  shelf  extends  20  kilometers  into  the  basin 
and  attains  a  maximum  thickness  of  1  kilometer.  In  the 
southern  part  of  the  basin,  the  thickness  increases  and  may 
exceed  3  kilometers. 

Next,  the  group  velocity  for  the  two  layer  model  was 
studied  for  variations  with  thickness.  The  resulting  curves 
are  shown  in  Figure  5.13.  The  near  source  velocities  used 
in  the  raytracing  may  be  attained  easily  by  thinning  the 
layer.  The  phase  velocities,  which  are  not  presented, 
actually  approach  a  constant  value  more  quickly,  as  the 
layer  thins,  than  the  group  velocities.  This  confirms  one 
of  the  assumptions  that  went  into  the  phase  velocity 
measurements:  that  the  near  source  phase  velocities  were 
constant  over  the  frequency  range  being  measured. 
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Figure  5.13  Comparison  of  theoretical  group  velocities  for  a  single  layer  over  a  half 
space  with  measured  values.  These  curves  illustrate  the  effect  of 
increasing  the  layer  thickness  (thicknesses  are  in  kilometers) . 
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Figure  5.13  illustrates  that  the  low  velocities 
modelled  from  raytracing  cannot  be  obtained  with  this  two 
dimensional  layered  model.  As  the  top  layer  is  thickened, 
the  group  velocity  curve  flattens,  but  no  velocities  lower 
than  approximately  2.4  km/s  occur.  If  the  structure 
continues  to  thicken,  the  group  velocity  will  eventually 
approach  a  constant  value  of  .9l94*V3j^,  or  2.52  km/s.  This 
is  because  the  depth  to  which  a  surface  wave  samples  the 
crust  is  on  the  order  of  1  wavelength.  Although  the 
wavelengths  are  not  known,  they  can  be  estimated  from 
average  velocities  and  frequency.  A  l  Hz  wave  with  a  phase 
velocity  of  2.5  km/s  has  a  wavelength  of  2.5  kilometers,  so 
as  the  basin  thickens  beyond  3  kilometers,  the  velocities 
will  be  controlled  almost  solely  by  the  Animikie  Group 
material  properties.  These  properties,  as  modelled  for  the 
average  Animikie  structure,  cannot  produce  the  low 
velocities  called  for  by  the  raytracing  model.  A  more 
complicated  model  is  required. 

In  the  final  raytracing  model,  most  of  the  rays  cross 
the  southern  basin  boundary  in  a  zone  where  the  velocities 
are  between  2.2  and  2.1  km/s.  By  changing  some  of  the 
properties  of  the  two  layers,  velocities  as  low  as  2.1  km/s 
could  be  generated  for  a  1  Hz  wave.  The  new  velocity  and 
density  structure  could  not  be  used,  however,  to  obtain  the 
observed  dispersion  curves.  The  parameters  for  the  low 
velocity  structure  are  shown  in  Figure  5.14.  The  shear  wave 
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5 . 14  Two  layer  model  which  yields  a  low  ^roup 
velocity  (2  km/s)  at  1  Hz 
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velocities  were  lowered,  especially  in  the  upper  layer. 

This  corresponds  to  increasing  Poisson's  ratio  or  decreasing 
the  rigidity  in  that  layer,  since  the  P-wave  velocity 
remains  unchanged.  Mooney  and  Bolt  (1965)  used  a  Poisson's 
ratio  of  0,33  when  modelling  a  sedimentary  basin  overlying 
basement.  Unfortunately,  this  is  difficult  to  justify 
geologically,  since  observed  deformation  increases  in  the 
part  of  the  basin  which  requires  low  velocities  (Southwick 
and  others,  1988) .  An  increase  in  deformation  implies 
increasing  metamorphic  grade  and  increasing  rigidity.  The 
density  was  also  lowered  in  the  upper  layer.  This  can  be 
justified  by  a  decrease  in  the  relative  amount  of  iron 
formation  present  as  the  total  amount  of  sedimentary 
material  increases. 

The  half  space  used  to  obtain  the  low  velocity  results 
has  a  slightly  higher  density,  which  could  correspond  to  a 
transition  in  the  material  underlying  the  basin  from  the 
Archean  greenstones  and  granites  to  either  the  rocks  of  the 
Great  Lakes  Tectonic  Zone,  or  the  Cuyuna  Range  rocks  of  the 
fold  and  thrust  terrane.  Figure  1.5  (Carlson,  1985).  The  S- 
wave  velocity  was  also  lowered  slightly,  which  is  more 
likely  to  correspond  to  the  Cuyuna  Range  rocks  since  they 
include  sedimentary  material. 

A  three  layer  model  might  be  used  to  match  the 
estimated  dispersion  curve  for  the  Animikie  basin,  however, 
adding  a  third  layer  introduces  4  additional  variables. 
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Also,  the  three  layer  model  cannot  be  readily  interpreted  in 
the  regions  where  the  Animikie  Group  has  been  tightly 
folded.  A  three  layer  model  for  the  Animikie  Group  might  be 
used  to  account  for  the  lateral  refraction  of  the  R1  arrival 
if  the  terrane  boundaries  were  redefined  and  the  refraction 
occurs  within  the  undeformed  portion  of  the  Animikie  basin. 
This  would  require  a  new  raytracing  model. 

The  difficulty  encountered  in  trying  to  reconcile  the 
raytracing  results  to  theoretical  dispersion  curves  makes 
the  case  for  the  regeneration  of  R1  within  the  fold  and 
thrust  terrane  stronger.  If  some  feature  within  that 
terrane  is  scattering  the  surface  waves,  it  is  not  necessary 
to  bend  the  raypaths  from  the  eastern  source  areas  to  get  R1 
to  arrive  at  the  array.  The  correlation  coefficients 
measured  at  Station  1  should  be  high  for  all  source  areas  if 
Rl  is  regenerated  and  this  is  not  the  case  for  Erie  and 
Reserve,  which  are  also  the  source  areas  which  would  be  most 
affected  by  the  low  velocity  structure.  Also,  if  the 
surface  waves  are  being  regenerated  within  the  fold  and 
thrust  terrane  by  interaction  with  a  specific  geologic 
feature,  the  maximum  cross^correlation  measured  for  Rl 
between  source  areas  should  be  high  for  all  of  the  array 
stations.  This  is  not  the  case. 
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5 . 5  Summary  of  Modelling  Results 

Two  dimensional  raytracing  provides  a  feasible 
mechanism  for  modelling  deviations  from  a  direct  raypath  for 
the  R1  surface  wave.  The  raypaths  for  Rl  are  sensitive  to 
both  the  geometry  and  velocities  of  the  major  tectonic 
terranes.  A  velocity  model  was  obtained  for  the  major 
Precambrian  tectonic  terranes  which  matches  the  observed 
group  travel  times  at  the  array  for  a  1  Hz  surface  wave 
generated  at  the  Mesabi  Range  mines.  The  velocity  model 
fits  previously  determined  P-wave  and  surface  wave 
velocities  and  accounts  for  the  bending  of  the  raypaths 
which  is  necessary  in  order  for  waves  generated  at  the 
easternmost  source  areas  to  reach  the  array. 

A  two  layer  structure  has  been  modelled  which  fits  the 
average  dispersion  curve  calculated  for  the  Animikie  basin. 
The  average  thickness  for  the  uppermost  layer  as  modelled  is 
1.2  kilometers.  This  two  layer  model  cannot  be  used  to 
account  for  the  low  velocities  required  by  the  velocity 
model  obtained  from  raytracing.  This  indicates  that  the 
southern  Animikie  basin  has  different  material  properties 
than  the  northern  part  of  the  basin  and  requires  a  more 
complicated  model  to  obtain  the  low  velocities  that  have 
been  implied  for  that  region.  It  may  also  be  possible  to 
match  the  velocities  by  using  a  three  or  four  layer  model, 
which  was  not  attempted. 
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5 . 6  Summary  and  Conclusions 

The  travel  times  measured  from  multiple  filter  analysis 
show  that  Rl  can  be  modelled  as  a  plane  wavefront  which 
propagates  from  northwest  to  southeast  across  the  array. 

The  wedge-like  structure  which  underlies  the  array  affects 
the  observed  group  velocities  and  dispersion.  Frequencies 
above  1.5  Hz  propagate  at  a  constant  group  velocity  of 
approximately  1.8  km/s.  Lower  frequencies  initially 
propagate  at  group  velocities  higher  than  1.8  km/s,  so 
dispersion  increases.  As  the  sedimentary  layer  thickens, 
the  low  frequencies  are  slowed  relative  to  high  frequencies 
and  a  subsequent  decrease  in  dispersion  occurs.  The  average 
group  velocity  as  a  function  of  frequency  for  all  six  array 
stations  can  be  modelled  as  a  sedimentary  layer,  2.1 
kilometers  thick,  overlying  a  half  space. 

The  azimuth  of  Rl  appears  to  be  solely  controlled  by 
the  low  velocity  material  which  underlies  the  array.  The 
results  from  raytracing  indicate  that  the  angle  at  which  Rl 
crosses  the  array  is  practically  independent  of  the  angle  of 
incidence  at  the  western  edge  of  the  Keweenawan  basin.  The 
azimuth  measured  from  phase  travel  times  had  a  constant 
value  of  320  degrees  for  all  frequencies  and  station 
combinations.  The  azimuth  measured  from  group  velocities 
averaged  320  degrees.  The  group  velocities  and  azimuths 
obtained  from  group  travel  times  were  higher  for  the 
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westernmost  stations  at  the  lowest  frequencies  (0.4  -  0.8 
Hz) .  This  indicates  the  presence  of  a  velocity  gradient  due 
to  the  wedge- like  structure  beneath  the  array. 

Travel  times  were  measured  to  Station  1  of  the  array 
for  Rl.  The  group  velocity  was  computed  assuming  a  straight 
raypath  to  Station  1.  The  velocities  were  very  consistent 
between  source  areas  and  ranged  from  2.65  km/s  to  2.85  km/s. 
The  consistency  and  narrow  range  of  measured  velocities 
suggest  that  the  true  raypath  for  Rl  is  close  to  a  straight 
raypath.  The  highest  group  velocities  occur  for  blasts  at 
U.S.  Steel.  Near  source  phase  velocity  measurements 
confirmed  the  presence  of  a  velocity  high  associated  with 
U.S.  Steel.  Although  the  straight  path  is  approximately 
correct,  some  lateral  refraction  must  occur,  presumably 
within  the  Animikie  basin,  in  order  for  the  surface  waves 
generated  at  Erie  and  Reserve  to  reach  the  array. 

Two  dimensional  raytracing  was  used  to  match  the  Rl 
travel  times  to  the  array.  A  velocity  model  for  the 
Animikie  basin  in  which  the  velocity  decreases  from  2.8  km/s 
at  the  source  areas  to  1.9  km/s  near  the  southeastern  margin 
of  the  basin  was  obtained.  The  raypaths  obtained  for  this 
structure  have  fairly  small  deviations  from  the  straight 
raypaths. 

An  average  group  dispersion  curve  for  the  Animikie 
basin  was  estimated  from  the  raytracing  model.  A  two  layer 
model  was  fit  to  the  curve.  A  thickness  for  the  Animikie 
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Group  of  1.2  kilometers  was  obtained.  To  obtain  velocities 
with  values  as  low  as  2.0  km/s,  the  shear  wave  velocity 
within  the  Animikie  Group  must  decrease  to  the  southeast. 

The  true  source  location  for  R1  could  not  be  determined 
unambiguously.  This  is  because  of  the  refraction  which 
occurs  as  Rl  arrives  at  the  array.  Rl  may  be  regenerated 
within  the  fold  and  thrust  terrane.  The  high  correlation 
coefficients  measured  at  Station  l,  and  results  from 
raytracing  suggest  this  is  a  possibility.  If  this  is  the 
case,  the  correlation  coefficients  should  be  high  at  all 
stations  of  the  array.  Since  the  correlation  coefficients 
are  not  high  between  records  from  different  sources  at 
Stations  2-6,  the  preferred  interp’*  -  -..ion  is  that  Rl 
propagates  as  a  wavefront  from  the  mines  to  the  array. 

The  results  for  Rl  indicate  that  mine  blasts  can  be 
very  useful  for  determining  regional  crustal  properties. 

The  location  of  the  mines  provided  good  coverage  of  the 
Animikie  basin.  The  correlation  between  sources  at  Station 
1  indicates  that  seismometer  location  is  very  important,  and 
that  array  processing  is  possible  if  the  instruments  are 
emplaced  over  basement. 

R2  was  modelled  as  a  phase  which  is  scattered  by  a 
point  source  near  the  Douglas  fault.  It  crosses  the  array 
at  a  velocity  of  1.2  km/s,  about  40  percent  slower  than  Rl. 
Since  anisotropy  and  wedge  structure  cannot  account  for  a 
velocity  difference  as  large  as  40  percent,  R2  must  be  a 
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different  mode.  R2  cannot  be  a  higher  Rayleigh  mode  unless 
there  is  a  velocity  inversion  within  the  basin,  since  higher 
modes  generally  travel  faster  than  the  fundamental  mode  at  a 
given  frequency.  R2  may  be  a  mixed  Love  and  Rayleigh  mode. 

The  identification  of  R2  is  difficult  for  Stations  1, 

2,  and  4.  R2  is  easily  identified  on  Stations  3,  5,  and  6. 
Modelling  of  R2  propagation  across  the  array  might  improve 
our  ability  to  pick  it  on  the  jther  array  stations.  The 
modal  nature  of  R2  and  the  velocity  structure  which  gives 
rise  to  its  scattering  might  be  determined  by  normal  mode 
synthesis  of  Lg  as  described  by  Maupin  (1989)  who  modelled 
Rg  and  Lg  propagation  across  the  North  Sea  Central  Graben. 


201 


References  Cited 


Aki,  K. ,  and  Richards,  P.G.,  1980,  Quantitative  Seismology: 
W.H.  Freeman  and  Company,  932  p. 

Bath,  M. ,  1974,  Spectral  Analysis  in  Geophysics;  Elsevier, 
Amsterdam,  5€3  p. 

Bloch,  S.,  and  Hales,  A.L.,  1968,  New  techniques  for  the 
determination  of  surface  wave  phase  velocities: 

Bulletin  of  the  Seismological  Society  of  America, 

V.  58,  p.  1021-1034. 

Bracewell,  R. ,  1965,  The  Fourier  Transform  and  Its 
Applications:  McGraw-Hill,  381  p. 

Carlson,  K.E.,  1985,  A  combined  analysis  of  gravity  and 
magnetic  anomalies  in  east-central  Minnesota: 

M.S.  thesis.  University  of  Minnesota,  138  p. 

Cerveny,  V.,  Molotkov,  I. A.,  and  Psencik,  I.,  1977,  Ray 
Method  in  Seismology,  Univerzita  Karlova,  Praha, 

214  p. 

Chandler,  V.W.,  1982,  Interpretation  of  high-resolution 
aeromagnetic  data  in  northeastern  Minnesota:  in 
Society  of  Exploration  Geophysicists,  52nd  Annual 
Meeting  and  Exposition,  Dallas  Texas,  Technical 
Program  Abstracts  and  Biographies,  p.  288-291. 

- ,  1983a,  Aeromagnetic  map  of  Minnesota,  Carlton  and  Pine 

Counties:  Minnesota  Geological  Survey, 

Aeromagnetic  Map  Series  Map  A-3,  scale  1:250,000. 

- ,  1983b,  Aeromagnetic  map  of  Minnesota,  east-central 

region:  Minnesota  Geological  Survey,  Aeromagnetic 
Map  Series  Map  A-4,  scale  1:250,00. 

- ,  1985,  Interpretation  of  Precambrian  geology  in 

Minnesota  using  low-altitude,  high  resolution 
aeromagnetic  data,  in  Hinze,  W.J.,  ed..  The 
Utility  of  Regional  Gravity  and  Magnetic  Anomaly 
Maps:  Society  of  Exploration  Geophysicists,  p. 

375-391. 

- ,  McSwiggen,  P.L.,  Morey,  G.B.,  Hinze,  W. J. ,  and 

Anderson,  R.R.,  1989,  Interpretation  of  seismic 
reflection,  gravity,  and  magnetic  data  across 
Middle  Proterozoic  Mid-Continent  rift  system, 
northwestern  Wisconsin,  eastern  Minnesota,  and 


central  Iowa:  American  Association  of  Petroleum 
Geologists  Bulletin,  v.  73,  p.  261-275. 


202 


Cline,  A.K.,  1974,  Six  subprograms  for  curve  fitting  using 
splines  under  tension:  Communications  of  the 
Association  for  Computing  Machinery,  V.  17,  p. 

220-223. 

Craddock,  C. ,  1972,  Keweenawan  geology  of  east-central  and 
southeastern  Minnesota,  in,  Sims,  P.K.,  and  Morey, 

G.B.,  eds..  Geology  of  Minnesota:  A  Centennial 
Volume:  Minnesota  Geological  Survey,  p.  416-424 

Der,  Z.A.,  Marshall,  M.E.,  O'Donnell,  A.,  and  McElfresh, 
T.W.,  1984,  Spatial  coherence  structure  and 
attenuation  of  the  Lg  phase,  site  effects,  and  the 
interpretation  of  the  Lg  coda:  Bulletin  of  the 
Seismological  Society  of  America,  v.  74,  p.  1125- 
1147. 

Dziewonski,  A.,  Bloch,  S.,  and  Landisman,  M. ,  1969,  A 
technique  for  the  analysis  of  transient  seismic 
signals:  Bulletin  of  the  Seismological  Society  of 

America,  v.  59,  p.  427-444. 

Dziewonski,  A.M.,  and  Hales,  A.L.,  1972,  Numerical  analysis 
of  dispersed  seismic  waves,  in  Bolt,  B.A. ,  ed. , 

Methods  in  Computational  Physics:  v.  ii. 

Seismology:  Surface  Waves  and  Earth  Oscillations, 
Academic  Press,  p.  39-86. 

Ewing,  W.M.,  Jardetzky,  W.S.,  and  Press,  F.,  1957,  Elastic 
Waves  in  Layered  Media:  McGraw-Hill,  380  p. 

Farnbach,  J.S.,  1975,  The  complex  envelope  in  seismic  signal 
analysis:  Bulletin  of  the  Seismological  Society 
of  America,  v.65,  p.  951-962. 

Farnham,  P.R.,  1967,  Crustal  structure  in  the  Keweenawan 
province  of  east-central  Minnesota  and  western 
Wisconsin:  Ph.D.  thesis.  University  of  Minnesota, 

464  p. 

Ferderer,  R.J.,  Jr.,  1988,  Werner  deconvolution  and  its 
application  to  the  Penokean  orogen,  east-central 
Minnesota:  Ph.D.  thesis.  University  of  Minnesota, 

284  p. 

Greenhalgh,  S.A.,  1979,  Studies  with  a  small  seismic  array 
in  east-central  Minnesota:  Ph.D.  thesis. 

University  of  Minnesota,  323  p. 


203 


Gregersen,  S.,  1984,  Lg-wave  propagation  and  crustal 
structure  differences  near  Denmark  and  the  North 
Sea:  Geophysical  Journal  of  the  Royal 
Astronomical  Society,  v.  79,  p.  217-234. 

Herrmann,  R. ,  ed.,  1978,  Computer  programs  in  earthquake 
seismology:  v.  II,  St.  Louis  University. 

Holst,  T.B.,  1982,  Evidence  for  multiple  deformation  during 
the  early  Proterozoic  Penokean  orogeny  in  the 
Middle  Precambrian  Thomson  Formation,  Minnesota: 
Canadian  Journal  of  Earth  Sciences,  v.  19,  p. 

2043-2047. 

- ,  1984,  Evidence  for  nappe  development  during  the  early 

Proterozoic  Penokean  orogeny,  Minnesota:  Geology, 

V.  12,  p.  135-138. 

Kennett,  B.L.N.,  1984,  Guided  wave  propagation  in  laterally 
varying  media-I.  theoretical  development: 

Geophysical  Journal  of  the  Royal  Astronomical 
Society,  v.  79,  p. 235-256. 

- ,  and  Mykkeltveit,  S.,  1984,  Guided  wave  propagation  in 

laterally  varying  media-II.  Lg-waves  in  north¬ 
western  Europe:  Geophysical  Journal  of  the  Royal 
Astronomical  Society,  v.79,  p.  257-269. 

King,  E.R. ,  and  Zietz,  I.,  1971,  Aeromagnetic  study  of  the 
Midcontinent  gravity  high  of  central  United 
States:  Bulletin  of  the  Geological  Society  of 
America,  v.  82,  p. 2187-2208. 

Knopoff,  L. ,  1969,  Phase  and  group  slowness  in  inhomogeneous 
media:  Journal  of  Geophysical  Research,  v.  74,  p. 

1701. 

- ,  Berry,  M.J.,  and  Schwab,  F.A.,  1967,  Tripartite  phase 

velocity  observations  in  laterally  heterogeneous 
regions:  Journal  of  Geophysical  Research,  v.  72, 
p.  2595-2602. 

Lee,  Y.W. ,  1960,  Statistical  Theory  of  Communication:  John 
Wiley  &  Sons,  Inc.,  509  p. 

Luetgart,  J.H.,  1988,  User's  manual  for  RAy84/R83PLT 
interactive  two-dimensional  raytracing/ synthetic 
seismogram  package:  U.S.  Geological  Survey  Open- 
File  Report  88-238,  52  p. 


204 


Maupin,  V.,  1989,  Numerical  modelling  of  Lg  wave  propagation 
across  the  North  Sea  Central  Graben:  Geophysical 
Journal  International,  v.  99,  p.  273-283. 

McSwiggen,  P.L.,  1987,  New  model  of  the  Midcontinent  rift  in 
eastern  Minnesota  and  western  Wisconsin: 

Tectonics,  v.  6,  p.  677-685. 

Mooney,  H.M.,  and  Bolt,  B.A.,  1965,  Dispersion  tables  for 
Rayleigh  waves  on  a  single  surface  layer:  VESIAC 
Special  Report  4410-102-X,  University  of  Michigan, 

198  p. 

- ,  Craddock,  C.,  Farnham,  P.R.,  Johnson,  S.H.,  and  Volz, 

G.,  1970,  Refraction  seismic  investigation  of  the 
northern  Midcontinent  gravity  high:  Journal  of 
Geophysical  Research,  v.  75,  p.  5056-5086. 

- ,  and  Walton,  1980,  Seismicity  and  tectonic 

relationships  for  upper  Great  Lakes  Precambrian 
shield  province:  Final  Report  NUREG/Cr-1569 ,  U.S. 
Nuclear  Regulatory  Commission,  85  p. 

Morey,  G.B.,  1972,  Petrology  of  Keweenawan  sandstones  in  the 
subsurface  of  southeastern  Minnesota,  in  Sims, 

P.K.,  and  Morey,  G.B.,  eds..  Geology  of  Minnesota: 

A  Centennial  Volume:  Minnesota  Geological  Survey, 
p.  436-449. 

- ,  1983,  Lower  Proterozoic  stratified  rocks  and  the 

Penokean  orogeny  in  east-central  Minnesota,  in 
Medaris,  L.G. ,  Jr.,  ed..  Early  Proterozoic  Geology 
of  the  Great  Lakes  Region:  Geological  Society  of 
America  Memoir  160,  p.  97-112. 

- ,  and  Mudrey,  M.G.,  Jr.,  1972,  Keweenawan  volcanic  rocks 

in  east-central  Minnesota,  in  Sims,  P.K.,  and 
Morey,  G.B.,  eds..  Geology  of  Minnesota:  A 
Centennial  Volume:  Minnesota  Geological  Survey, 
p.  425-430. 

- ,  and  Ojakangas,  R.W.,  1982,  Keweenawan  sedimentary 

rocks  of  eastern  Minnesota  and  northwestern 
Wisconsin,  in  Wold,  R. J. ,  and  Hinze,  W.J.,  eds.. 

Geology  and  Tectonics  of  the  Lake  Superior  Basin: 
Geological  Society  of  America  Memoir  160,  p.  135- 
146. 

- ,  Sims,  P.K.,  Cannon,  W.F.,  Mudrey,  M.G.,  Jr.,  and 

Southwick,  D.L.,  1982,  Geologic  map  of  the  Lake 
Superior  region,  bedrock  geology:  Minnesota 


Geological  Survey  State  Map  Series  S-13,  scale 

1:1,000,000. 


205 


Mosher,  C.C.,  1980,  Signal  processing  techniques  applied  to 
a  small  circular  seismic  array:  Ph.D.  thesis. 
University  of  Minnesota,  267  p. 

Nafe,  J.E.,  and  Brune,  J.N.,  1960,  Observations  of  phase 
velocity  for  Rayleigh  waves  in  the  period  range 
100  to  400  seconds:  Bulletin  of  the  Seismological 
Society  of  America,  v.  50,  p.  427-440. 

Neidel,  N.S.,  and  Taner,  M.T.,  1971,  Semblance  and  other 
coherency  measures  for  multichannel  data: 

Geophysics,  v.  36,  p.  482-497. 

Oliver,  J.,  Ewing,  M.,  and  Press,  F.,  1955,  Crustal 

structure  of  the  Arctic  regions  from  the  Lg  phase: 
Bulletin  of  the  Geological  Society  of  America,  v. 

66,  p.  1063-1074. 

Otnes,  R.K.,  and  Enochsen,  L. ,  1978,  Applied  Time  Series 
Analysis:  John  Wiley  and  Sons,  449  p. 

Press,  F.,  and  Ewing,  M. ,  1952,  Two  slow  surface  waves 
across  North  America:  Bulletin  of  the 
Seismological  Society  of  America,  v.  42,  p.  219-228. 

Richter,  C.F.,  1958,  Elementary  Seismology:  W.H.  Freeman 
and  Company,  768  p. 

Southwick,  D.L.,  Meyer,  G.N.,  and  Mills,  S.J.,  1986, 
Scientific  core  drilling  in  central  Minnesota: 

Summary  of  lithologic  and  geochemical  results: 

Minnesota  Geological  Survey  Information  Circular 
23,  186  p. 

- ,  Morey,  G.B.,  and  McSwiggen,  P.L.,  1988,  Geologic  map 

(scale  1:250,000)  of  the  Penokean  orogen,  central 
and  eastern  Minnesota,  and  accompanying  text: 

Minnesota  Geological  Survey  Report  of 
Investigations  37,  25  p. 

Wang,  C.Y.,  1985a,  Computer  Programs  in  Seismology 
Intro(III):  St.  Louis  University,  14  p. 

- ,  1985b,  Computer  Programs  in  Seismology  Surface85 (III) : 

St.  Louis  University,  5  p. 

— -,  1985c,  Computer  Programs  in  Seismology  Reign85 (III)  : 
St.  Louis  University,  4  p. 


206 


- ,  1985d,  Computer  Programs  in  Seismology  DpegnSS (III) : 

St.  Louis  University,  2  p. 

Wright,  H.E.,  Mattson,  L.A. ,  and  Thomas,  j.,  1970,  Geology 
of  the  Cloquet  quadrangle,  Carlton  county, 

Minnesota;  Minnesota  Geological  Survey  General 
Map  Series  Map  GM-3. 


D1 

■H 


4 


